Commit dca2dca0 authored by Adam Paquette's avatar Adam Paquette
Browse files

Made initial additions of voronoi to the code base.

parent 6ebb086c

autocnet/cg/cg.py

+69 −71
Changes for autocnet/cg/cg.py: 69 added lines, 71 removed lines.
Original line number Diff line number Diff line 
import pandas as pd

import numpy as np

import geopandas as gpd

import ogr



from scipy.spatial import ConvexHull

from scipy.spatial import Voronoi

from shapely.geometry import Polygon, Point

from shapely.affinity import scale

import cv2



from autocnet.utils import utils
@@ -150,102 +153,94 @@ def vor(graph, clean_keys=[], s=30): 
            if neighbors != num_neighbors:

                raise AssertionError('The graph is not complete')



        source = graph.edges(0)[0]

        destination = graph.edges()[0][1]

        edge = graph.edge[source][destination]

        matches, _ = edge.clean(clean_keys=['fundamental'])

        kps = edge.source.get_keypoint_coordinates(index=matches['source_idx'], homogeneous=True)

        intersection, proj_gdf, source_gdf = compute_intersection(graph.nodes()[0], graph, clean_keys)



        intersection_poly = compute_intersection(graph, clean_keys)

        intersection_points = intersection_poly.GetGeometryRef(0).GetPoints()

        intersection_ind = intersection.query('overlaps_all == True').index.values[0]



        centroid = intersection_poly.Centroid().GetPoint()

        source_node = graph.nodes()[0]

        for e in graph.edges():

            if e[0] != source_node:

                source_node = e[0]

                intersection, proj_gdf, source_gdf = compute_intersection(source_node, graph, clean_keys)

                intersection_ind = intersection.query('overlaps_all == True').index.values[0]



        points = np.asarray(kps)

        voronoi_np = []

            edge = graph.edge[e[0]][e[1]]

            matches, mask = edge.clean(clean_keys = clean_keys)

            kps = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)



        point_cloud = ogr.Geometry(ogr.wkbMultiPoint)

        for p in points:

            point = ogr.Geometry(ogr.wkbPoint)

            point.AddPoint(double(p[0]), double(p[1]))

            point_cloud.AddGeometry(point)

            kps['geometry'] = kps.apply(lambda x: Point(x['x'], x['y']), axis = 1)

            kps.mask = kps['geometry'].apply(lambda x: intersection.geometry.contains(x).all())



        intersection_cloud = intersection_poly.Intersection(point_cloud)

            # Creates a mask for displaying the voronoi points

            # Currently erronious and produces NaN values

            # mask[mask] = kps.mask

            # edge.masks = ('voronoi', mask)



        for p in intersection_cloud:

            point = p.GetPoint(0)

            voronoi_np.append(point)

            keypoints = []

            kps[kps.mask].apply(lambda x: keypoints.append((x['x'], x['y'])), axis = 1)

            coords = source_gdf.geometry.apply(lambda x:scale(x, s, s).exterior.coords)



        keypoints = np.asarray(voronoi_np)

        voronoi_pd = pd.DataFrame(data=voronoi_np, columns=['x', 'y', 'homogenious'])

            for point in coords:

                keypoints = np.vstack((keypoints, point))



        # Based on the keypoints found in the method

        inters = np.empty((len(intersection_points),2))



        for g, (i, j) in enumerate(intersection_points):

            scaledx, scaledy = utils.scale_point((i, j), centroid, s)

            point = np.array([scaledx, scaledy])

            inters[g] = point



        keypoints = np.vstack((keypoints[:,:2], inters))

            vor = Voronoi(keypoints)

            voronoi_df = pd.DataFrame(data = kps, columns = ['x', 'y', 'weight'])



        voronoi_df = pd.DataFrame(columns=["x", "y", "weights"])

        voronoi_df["x"] = kps['x']

        voronoi_df["y"] = kps['y']



        poly_array = []

        for i, region in enumerate(vor.regions):

            region_point = vor.points[np.argwhere(vor.point_region==i)]

            i = 0

            vor_points = np.asarray(vor.points)

            for region in vor.regions:

                region_point = vor_points[np.argwhere(vor.point_region==i)]

                if -1 not in region:

                    polygon_points = [vor.vertices[i] for i in region]

                    if len(polygon_points) != 0:

                    polygon = utils.array_to_poly(polygon_points)

                    intersection = polygon.Intersection(intersection_poly)

                    poly_array = np.append(poly_array, intersection)

                    polygon_area = intersection.GetArea()

                        polygon = Polygon(polygon_points)

                        poly_area = polygon.intersection(intersection.geometry[intersection_ind]).area

                        voronoi_df.loc[(voronoi_df["x"] == region_point[0][0][0]) &

                                       (voronoi_df["y"] == region_point[0][0][1]),

                                   'vor_weights'] = polygon_area

                                       'weight'] = poly_area

                i += 1



            edge.weight['vor_weight'] = voronoi_df['weight']



        return vor, voronoi_df



def compute_intersection(source, graph, clean_keys=[]):

    if type(source) is int:

        source = graph.node[source]



def compute_intersection(graph, clean_keys=[]):

    source_num = graph.edges()[0][0]

    source = graph.node[source_num]

    source_corners = source.geodata.xy_corners

    total_intersect_poly = utils.array_to_poly(source_corners)

    orig_poly = utils.array_to_poly(source_corners)

    source_poly = Polygon(source_corners)



    source_gdf = gpd.GeoDataFrame({'source_geom': [source_poly], 'source_node': [source.node_id]}).set_geometry('source_geom')



    proj_list = []

    proj_nodes = []

    # Begin iterating through the nodes in the graph excluding the source node

    for n in graph.nodes_iter():

        if n == source_num:

        if n == source.node_id:

            continue



        # Define the edge, matches, and destination based on the zero node and the nth node

        edge = graph.edge[source_num][n]

        matches, _ = edge.clean(clean_keys=clean_keys)

        destination = edge.destination

        destination = graph.node[n]

        destination_corners = destination.geodata.xy_corners



        kp1 = edge.source.get_keypoint_coordinates(index=matches['source_idx'], homogeneous=True)

        kp2 = edge.destination.get_keypoint_coordinates(index=matches['destination_idx'], homogeneous=True)

        edge = graph.edge[source.node_id][destination.node_id]



        # If the source image has coordinate transformation data, us the footprint and

        # coordinate transforms as it will produce a more accurate intersection

        if edge.source.geodata.coordinate_transformation.this is not None:

            source_footprint_poly = edge.source.geodata.footprint

            destination_footprint_poly = edge.destination.geodata.footprint

        # Will still need the check but the helper function will make these calls much easier to understand

        if source.node_id > destination.node_id:

            kp2 = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)

            kp1 = edge.get_keypoints('destination', clean_keys = clean_keys, homogeneous = True)

        else:

            kp2 = edge.get_keypoints('destination', clean_keys = clean_keys, homogeneous = True)

            kp1 = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)



            intersection_poly = destination_footprint_poly.Intersection(source_footprint_poly)

        # If the source image has coordinate transformation data

        if (source.geodata.coordinate_transformation.this != None) and \

        (destination.geodata.coordinate_transformation.this != None):

            proj_poly = Polygon(destination.geodata.xy_corners)



        # Else, use the homography transform to get an intersection of the two images

        else:

            H, mask = cv2.findHomography(kp2.values,

                                         kp1.values,

                                         cv2.RANSAC,

                                         2.0)



            H, mask = cv2.findHomography(kp2.values, kp1.values, cv2.RANSAC, 2.0)

            proj_corners = []

            for c in destination_corners:

                x, y, h = utils.reproj_point(H, c)
@@ -254,12 +249,15 @@ def compute_intersection(graph, clean_keys=[]): 
                h /= h

                proj_corners.append((x, y))



            proj_poly = utils.array_to_poly(proj_corners)

            proj_poly = Polygon(proj_corners)



            intersection_poly = orig_poly.Intersection(proj_poly)

        proj_list.append(proj_poly)

        proj_nodes.append(n)



        # Intersect the newly calculated intersection with the current total intersection

        # to get the new bound of the overlap area

        total_intersect_poly = intersection_poly.Intersection(total_intersect_poly)

    proj_gdf = gpd.GeoDataFrame({'proj_geom': proj_list, 'proj_node': proj_nodes}).set_geometry('proj_geom')



        return total_intersect_poly
 
    intersect_gdf = gpd.overlay(source_gdf, proj_gdf, how='intersection')

    intersect_gdf = intersect_gdf.rename(columns = {'geometry':'intersect_geom'}).set_geometry('intersect_geom')

    intersect_gdf['overlaps_all'] = intersect_gdf.geometry.apply(lambda x:proj_gdf.geometry.contains(scale(x, .9, .9)).all())

    intersection = intersect_gdf.query('overlaps_all == True').set_geometry('intersect_geom')

    return intersection, proj_gdf, source_gdf

autocnet/graph/edge.py

+3 −3
Changes for autocnet/graph/edge.py: 3 added lines, 3 removed lines.
Original line number Diff line number Diff line
@@ -465,11 +465,11 @@ class Edge(dict, MutableMapping): 
            node = node.lower()



        if type(node) is Node:

            node = node.node_id

            node = node['node_id']



        if node == "source" or node == "s" or node == self.source.node_id:

        if node == "source" or node == "s" or node == self.source['node_id']:

            return self.source.get_keypoint_coordinates(index=matches['source_idx'], **kwargs)

        if node == "destination" or node == "d" or node == self.destination.node_id:

        if node == "destination" or node == "d" or node == self.destination['node_id']:

            return self.destination.get_keypoint_coordinates(index=matches['destination_idx'], **kwargs)



    def decompose(self, maxiterations=3):

autocnet/graph/network.py

+3 −0
Changes for autocnet/graph/network.py: 3 added lines, 0 removed lines.
Original line number Diff line number Diff line
@@ -683,3 +683,6 @@ class CandidateGraph(nx.Graph): 
        """

        edges = [(u, v) for u, v, edge in self.edges_iter(data=True) if func(edge, *args, **kwargs)]

        return self.create_edge_subgraph(edges)



    def compute_vor_weight(self, clean_keys = [], **kwargs):

        pass

environment.yml

+2 −0
Changes for environment.yml: 2 added lines, 0 removed lines.
Original line number Diff line number Diff line
@@ -17,6 +17,8 @@ dependencies: 
  - numexpr

  - numpy>1.10.0

  - pandas

  - geopandas

  - shapely

  - pyyaml

  - scipy>0.17.0

  - sqlalchemy

notebooks/Testing_geopandas.ipynb

+54 −37
Changes for notebooks/Testing_geopandas.ipynb: 54 added lines, 37 removed lines.
Original line number Diff line number Diff line 
%% Cell type:code id: tags:



``` python

import os

import sys



sys.path.insert(0, os.path.abspath('..'))



import cv2

import numpy as np

import geopandas as gpd

import pandas as pd

import networkx as nx

from shapely.geometry import Polygon as Poly, Point

from shapely.affinity import scale

import ogr

import time

from scipy.spatial import Voronoi



from autocnet.examples import get_path

from autocnet.graph.network import CandidateGraph

from autocnet.cg import cg



from unittest.mock import Mock

from unittest.mock import MagicMock

from autocnet.graph import edge

from autocnet.graph import node

from autocnet.graph.node import Node

from plio.io import io_gdal



from autocnet.transformation.transformations import Homography

from autocnet.utils.utils import array_to_poly



from IPython.display import display



%pylab inline

figsize(20, 20)

```



%% Output



    Populating the interactive namespace from numpy and matplotlib



%% Cell type:code id: tags:



``` python

def scale_point(point, centroid, scalar):

    point = np.asarray(point)

    centroid = centroid[:2]

    vector = ((point - centroid)*scalar) + centroid

    return (vector)

```



%% Cell type:code id: tags:



``` python

def reproj_point(H, point):

    """

    Reproject a pixel in one image into another image



    Parameters

    ----------

    H : object

        (3,3) ndarray or Homography object



    corner : iterable

             A 2 element iterable in the form x, y

    """

    if len(point) == 2:

        coords = np.array([point[0],point[1],1])

    elif len(point) == 3:

        coords = np.asarray(point)

        coords *= coords[-1]  # Homogenize



    return H.dot(coords)

```



%% Cell type:code id: tags:



``` python

def compute_intersection(source, graph, clean_keys=[]):

    if type(source) is int:

        source = graph.node[source]



    source_corners = source.geodata.xy_corners

    source_poly = Poly(source_corners)



    source_gdf = gpd.GeoDataFrame({'source_geom': [source_poly], 'source_node': [source.node_id]}).set_geometry('source_geom')

    source_gdf = gpd.GeoDataFrame({'source_geom': [source_poly], 'source_node': [source['node_id']]}).set_geometry('source_geom')



    proj_list = []

    proj_nodes = []

    # Begin iterating through the nodes in the graph excluding the source node

    for n in graph.nodes_iter():

        if n == source.node_id:

        if n == source['node_id']:

            continue



        # Define the edge, matches, and destination based on the zero node and the nth node

        destination = graph.node[n]

        destination_corners = destination.geodata.xy_corners

        print(source, destination)



        edge = graph.edge[source.node_id][destination.node_id]

        edge = graph.edge[source['node_id']][destination['node_id']]



        # Will still need the check but the helper function will make these calls much easier to understand

        if source.node_id > destination.node_id:

        if source['node_id'] > destination['node_id']:

            kp2 = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)

            kp1 = edge.get_keypoints('destination', clean_keys = clean_keys, homogeneous = True)

        else:

            kp2 = edge.get_keypoints('destination', clean_keys = clean_keys, homogeneous = True)

            kp1 = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)



        # If the source image has coordinate transformation data

        if (source.geodata.coordinate_transformation.this != None) and \

        (destination.geodata.coordinate_transformation.this != None):

            proj_poly = Poly(destination.geodata.xy_corners)



        # Else, use the homography transform to get an intersection of the two images

        else:

            H, mask = cv2.findHomography(kp2.values, kp1.values, cv2.RANSAC, 2.0)

            proj_corners = []

            for c in destination_corners:

                x, y, h = reproj_point(H, c)

                x /= h

                y /= h

                h /= h

                proj_corners.append((x, y))



            proj_poly = Poly(proj_corners)



        proj_list.append(proj_poly)

        proj_nodes.append(n)



    proj_gdf = gpd.GeoDataFrame({'proj_geom': proj_list, 'proj_node': proj_nodes}).set_geometry('proj_geom')



    intersect_gdf = gpd.overlay(source_gdf, proj_gdf, how='intersection')

    intersect_gdf = intersect_gdf.rename(columns = {'geometry':'intersect_geom'}).set_geometry('intersect_geom')

    intersect_gdf['overlaps_all'] = intersect_gdf.geometry.apply(lambda x:proj_gdf.geometry.contains(scale(x, .9, .9)).all())

    intersection = intersect_gdf.query('overlaps_all == True').set_geometry('intersect_geom')

    return intersection, proj_gdf, source_gdf

```



%% Cell type:code id: tags:



``` python

def vor(graph, clean_keys=[], verbose=False, s = 30):



    num_neighbors = len(graph.nodes()) - 1

    for n in graph.nodes():

        neighbors = len(graph.neighbors(n))

        if neighbors != num_neighbors:

            raise AssertionError('The graph is not complete')



    # Have compute_intersection get the intersections of the subgraph, the projected nodes in the

    # subgraph, and the source image all in there own geopandas dataframes

    intersection, proj_gdf, source_gdf = compute_intersection(graph.nodes()[0], graph, clean_keys)



    # Get the intersection for the given subgrpah

    intersection_ind = intersection.query('overlaps_all == True').index.values[0]



    source_node = graph.nodes()[0]

    # Set the source, destination, edge, matches, and keypoints(kps) based on the first edge in the subgraph

    # which should be the smallest node in the graph

    for e in graph.edges():

        # Recompute the intersection if the source node of the n + 1 edge is different from the n edge

        if e[0] != source_node:

            source_node = e[0]

            intersection, proj_gdf, source_gdf = compute_intersection(source_node, graph, clean_keys)

            intersection_ind = intersection.query('overlaps_all == True').index.values[0]





        edge = graph.edge[e[0]][e[1]]

        matches, mask = edge.clean(clean_keys = clean_keys)

        kps = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)



        # Get all of the keypoints that are inside of the intersection then mask out the points

        # in kps that lie outside of the intersection

        # All of these operations on kps are extremely expensive

        kps['geometry'] = kps.apply(lambda x: Point(x['x'], x['y']), axis = 1)

        kps.mask = kps['geometry'].apply(lambda x: intersection.geometry.contains(x).all())



        # Creates a mask for displaying the voronoi points

        # Currently erronious and produces NaN values

        # mask[mask] = kps.mask

        # edge.masks = ('voronoi', mask)



        keypoints = []

        kps[kps.mask].apply(lambda x: keypoints.append((x['x'], x['y'])), axis = 1)

        coords = source_gdf.geometry.apply(lambda x:scale(x, s, s).exterior.coords)



        for point in coords:

            keypoints = np.vstack((keypoints, point))



        vor = Voronoi(keypoints)



        # For weight computation

        # Should move to its own method

        voronoi_df = pd.DataFrame(data = kps, columns = ['x', 'y', 'weight'])



        i = 0

        vor_points = np.asarray(vor.points)

        for region in vor.regions:

            region_point = vor_points[np.argwhere(vor.point_region==i)]

            if not -1 in region:

                polygon_points = [vor.vertices[i] for i in region]

                if len(polygon_points) != 0:

                    polygon = Poly(polygon_points)

                    poly_area = polygon.intersection(intersection.geometry[intersection_ind]).area

                    voronoi_df.loc[(voronoi_df["x"] == region_point[0][0][0]) &

                                   (voronoi_df["y"] == region_point[0][0][1]),

                                   'weight'] = poly_area

            i+=1



        edge.weight['vor_weight'] = voronoi_df['weight']

        edge['weights']['vor_weight'] = voronoi_df['weight']

```



%% Cell type:code id: tags:



``` python

#Point to the adjacency Graph

adjacency = get_path('vor_adjacency1.json')

basepath = get_path('Apollo15')

cang = CandidateGraph.from_adjacency(adjacency, basepath=basepath)



#Apply SIFT to extract f"eatures

cang.extract_features(method='sift', extractor_parameters={'nfeatures':1500})



#Match

cang.match_features()

cang.match()



#Apply outlier detection

# cg.apply_func_to_edges(Edge.symmetry_check)

cang.symmetry_checks()

# cg.apply_func_to_edges(Edge.ratio_check)

cang.ratio_checks()





#Compute a homography and apply RANSAC

cang.apply_func_to_edges("compute_fundamental_matrix", clean_keys=['ratio', 'symmetry'])

```



%% Output



    /home/acpaquette/autocnet/autocnet/transformation/fundamental_matrix.py:252: UserWarning: Unable to apply MLE.  Not enough correspondences.  Returning with a RANSAC computed F matrix.

      warnings.warn("Unable to apply MLE.  Not enough correspondences.  Returning with a RANSAC computed F matrix.")



%% Cell type:code id: tags:



``` python

def compute_weights(cang, source_id, clean_keys = [], verbose = False):

    cliques = compute_cliques(cang, source_id)

    source = cang.node[source_id]

    for g in cliques:

        subgraph = cang.create_node_subgraph(g)

        gdf_tup = vor(subgraph, source, clean_keys, verbose)

```



%% Cell type:code id: tags:



``` python

def compute_cliques(graph, source_id):

    if source_id is not None:

        return list(nx.cliques_containing_node(cang, nodes=source_id))

    else:

        return list(nx.find_cliques(graph))

```



%% Cell type:code id: tags:



``` python

start = time.clock()

cliques = compute_cliques(cang, None)

for g in cliques:

    subgraph = cang.create_node_subgraph(g)

    vor(subgraph, clean_keys = ['fundamental'], verbose = False)

# subgraph = cang.create_node_subgraph([1, 4])

# voronoi = vor(cang, clean_keys = ['fundamental'], verbose = False)

    vor(subgraph, clean_keys = ['fundamental'])

print(time.clock() - start)



start = time.clock()

for e in cang.edges_iter():

    subgraph = cang.create_node_subgraph(e)

    vor(subgraph, clean_keys = ['fundamental'])

print(time.clock() - start)

```



%% Output



    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    Getting CT

    CT <osgeo.osr.CoordinateTransformation; proxy of None >

    2.850580000000001

    2.9477230000000016



%% Cell type:code id: tags:



``` python

source_keypoint_df = pd.DataFrame({'x': (15, 18, 18, 12, 12), 'y': (6, 10, 15, 15, 10)})

destination_keypoint_df = pd.DataFrame({'x': (5, 8, 8, 2, 2), 'y': (1, 5, 10, 10, 5)})

keypoint_matches = [[0, 0, 1, 0],

                    [0, 1, 1, 1],

                    [0, 2, 1, 2],

                    [0, 3, 1, 3],

                    [0, 4, 1, 4]]



matches_df = pd.DataFrame(data = keypoint_matches, columns = ['source_image', 'source_idx',

                                                                  'destination_image', 'destination_idx'])

# Source and Destination Node Setup

source_node = MagicMock(spec = Node())

destination_node = MagicMock(spec = Node())



source_node.get_keypoint_coordinates = MagicMock(return_value=source_keypoint_df)

destination_node.get_keypoint_coordinates = MagicMock(return_value=destination_keypoint_df)

source_node.node_id = 0

destination_node.node_id = 1

source_node.__get_item__ = 0

destination_node.__get_item__ = 1



source_geodata = Mock(spec = io_gdal.GeoDataset)

destination_geodata = Mock(spec = io_gdal.GeoDataset)



source_node.geodata = source_geodata

destination_node.geodata = destination_geodata



source_node.geodata.coordinate_transformation.this = None

destination_node.geodata.coordinate_transformation.this = None



source_corners = [(0, 0),

                  (20, 0),

                  (20, 20),

                  (0, 20)]



destination_corners = [(0, 0),

                       (20, 0),

                       (20, 20),

                       (0, 20)]



source_node.geodata.xy_corners = source_corners

destination_node.geodata.xy_corners = destination_corners



# Edge Setup

e = edge.Edge(source = source_node, destination = destination_node)



e.clean = MagicMock(return_value=(matches_df, None))

e.matches = matches_df



def side_effect(node, clean_keys, **kwargs):

    if type(node) is str:

        node = node.lower()



    if isinstance(node, Node):

        node = node.node_id

        node = node['node_id']



    if node == "source" or node == "s" or node == source_node.node_id:

    if node == "source" or node == "s" or node == source_node['node_id']:

        return e.source.get_keypoint_coordinates().copy(deep=True)

    if node == "destination" or node == "d" or node == destination_node.node_id:

    if node == "destination" or node == "d" or node == destination_node['node_id']:

        return e.destination.get_keypoint_coordinates().copy(deep=True)



e.get_keypoints = MagicMock(side_effect = side_effect)



cang = MagicMock(spec = CandidateGraph())



cang.nodes = MagicMock(return_value=([0, 1]))

cang.node = [source_node, destination_node]

cang.nodes_iter = MagicMock(return_value=([0, 1]))

cang.neighbors = MagicMock(return_value = [1])

cang.edges = MagicMock(return_value=([(0, 1)]))

cang.edge = {0:{1:e}, 1:{0:e}}

```



%% Cell type:code id: tags:



``` python

vor(cang, clean_keys = ['fundamental'])

```



%% Output



    0    22.50

    1    26.25

    2    37.50

    3    37.50

    4    26.25

    Name: weight, dtype: float64

    <MagicMock spec='Node' id='140434391257040'> <MagicMock spec='Node' id='140434391257040'>



    ---------------------------------------------------------------------------

    KeyError                                  Traceback (most recent call last)

    <ipython-input-22-60125f4ee637> in <module>()

    ----> 1 vor(cang, clean_keys = ['fundamental'])



    <ipython-input-5-733ea7cc66d3> in vor(graph, clean_keys, verbose, s)

          9     # Have compute_intersection get the intersections of the subgraph, the projected nodes in the

         10     # subgraph, and the source image all in there own geopandas dataframes

    ---> 11     intersection, proj_gdf, source_gdf = compute_intersection(graph.nodes()[0], graph, clean_keys)

         12

         13     # Get the intersection for the given subgrpah

    <ipython-input-12-7fbccd493ab1> in compute_intersection(source, graph, clean_keys)

         20         print(source, destination)

         21

    ---> 22         edge = graph.edge[source['node_id']][destination['node_id']]

         23

         24         # Will still need the check but the helper function will make these calls much easier to understand

    KeyError: <MagicMock name='mock.__getitem__()' id='140434391038104'>



%% Cell type:code id: tags:



``` python

# adjacency = get_path('two_image_adjacency.json')

# basepath = get_path('Apollo15')

# cg = CandidateGraph.from_adjacency(adjacency, basepath=basepath)

keypoint_df = pd.DataFrame({'x': (15, 18, 18, 12, 12), 'y': (6, 10, 15, 15, 10)})

keypoint_matches = [[0, 0, 1, 0],

                    [0, 1, 1, 1],

                    [0, 2, 1, 2],

                    [0, 3, 1, 3],

                    [0, 4, 1, 4]]



matches_df = pd.DataFrame(data = keypoint_matches, columns = ['source_image', 'source_idx',

                                                                  'destination_image', 'destination_idx'])

e = edge.Edge(source = source_node, destination = destination_node)



e.clean = MagicMock(return_value=(matches_df, None))

e.matches = matches_df



source_node = MagicMock(spec = node.Node())

destination_node = MagicMock(spec = node.Node())



cang = MagicMock(spec = CandidateGraph())



cang.nodes = MagicMock(return_value=([0, 1]))

cang.node = [source_node, destination_node]

cang.nodes_iter = MagicMock(return_value=([0, 1]))

cang.neighbors = MagicMock(return_value=([0]))

cang.edges = MagicMock(return_value=([(0, 1)]))

cang.edge = {0:{1:e}, 1:{0:e}}



source_node.get_keypoint_coordinates = MagicMock(return_value=keypoint_df)

destination_node.get_keypoint_coordinates = MagicMock(return_value=keypoint_df)

source_node.node_id = 0

destination_node.node_id = 1

source_node['node_id'] = 0

destination_node['node_id'] = 1



def side_effect(node, clean_keys, **kwargs):

    if type(node) is str:

        node = node.lower()



    if isinstance(node, Node):

        node = node.node_id

        node = node['node_id']



    if node == "source" or node == "s" or node == source_node.node_id:

    if node == "source" or node == "s" or node == source_node['node_id']:

        return e.source.get_keypoint_coordinates()

    if node == "destination" or node == "d" or node == destination_node.node_id:

    if node == "destination" or node == "d" or node == destination_node['node_id']:

        return e.destination.get_keypoint_coordinates()



e.get_keypoints = MagicMock(side_effect = side_effect)



e.source = source_node

e.destination = destination_node



source_geodata = Mock(spec = io_gdal.GeoDataset)

destination_geodata = Mock(spec = io_gdal.GeoDataset)



e.source.geodata = source_geodata

e.destination.geodata = destination_geodata



source_corners = [(0, 0),

                  (20, 0),

                  (20, 20),

                  (0, 20)]



destination_corners = [(10, 5),

                       (30, 5),

                       (30, 25),

                       (10, 25)]



source_xy_extent = [(0, 20), (0, 20)]



destination_xy_extent = [(10, 30), (5, 25)]



source_poly = array_to_poly(source_corners)

destination_poly = array_to_poly(destination_corners)



vals = {(10, 5):(10, 5), (20, 5):(20, 5), (20, 20):(20, 20), (10, 20):(10, 20)}



def latlon_to_pixel(i, j):

    return vals[(i, j)]



e.source.geodata.latlon_to_pixel = MagicMock(side_effect = latlon_to_pixel)

e.destination.geodata.latlon_to_pixel = MagicMock(side_effect = latlon_to_pixel)



e.source.geodata.footprint = source_poly

e.source.geodata.xy_corners = source_corners

e.source.geodata.xy_extent = source_xy_extent

e.destination.geodata.footprint = destination_poly

e.destination.geodata.xy_corners = destination_corners

e.destination.geodata.xy_extent = destination_xy_extent

```



%% Cell type:code id: tags:



``` python

vor(cang, clean_keys = ['fundamental'])

```



%% Output



    0    22.50

    1    26.25

    2    37.50

    3    37.50

    4    26.25

    Name: weight, dtype: float64



%% Cell type:code id: tags:



``` python

```



%% Cell type:code id: tags:



``` python

def vor_plot(graph, clean_keys = []):



    num_neighbors = len(graph.nodes()) - 1

    for n in graph.nodes():

        neighbors = len(graph.neighbors(n))

        if neighbors != num_neighbors:

            raise AssertionError('The graph is not complete')



    intersection, proj_gdf, source_gdf = compute_intersection(graph.nodes()[0], graph, clean_keys)



    source_node = graph.nodes()[0]



    for e in graph.edges():

        if e[0] != source_node:

            source_node = e[0]

            start_inter = time.clock()

            intersection, proj_gdf, source_gdf = compute_intersection(source_node, graph, clean_keys)



        edge = graph.edge[e[0]][e[1]]

        kps = edge.get_keypoints('source', clean_keys = clean_keys, homogeneous = True)

        kps = kps[edge.vor_mask]



        vor = edge.voronoi



        i = 0

        poly_array = []

        vor_points = np.asarray(vor.points)

        for region in vor.regions:

            region_point = vor_points[np.argwhere(vor.point_region==i)]

            if not -1 in region:

                polygon_points = [vor.vertices[i] for i in region]

                if len(polygon_points) != 0:

                    polygon = Poly(polygon_points)

                    poly_array.append(polygon)

            i+=1



        poly_gdf = gpd.GeoDataFrame(data = poly_array, columns = ['geometry'])

        vor_poly_gdf = gpd.overlay(poly_gdf, intersection, how='intersection')



        ax = proj_gdf.query('proj_node != ' + str(n)).plot(color='#388aff', alpha=0.1)

        proj_gdf.query('proj_node == ' + str(n)).plot(color = 'g', alpha=0.1, ax = ax)

        source_gdf.plot(color='r', alpha=0.1, ax=ax)

        vor_poly_gdf.plot(cmap='Set1', alpha=.5, ax=ax)

#             graph.edge[0][3].plot(clean_keys = ['fundamental'], ax = ax)

        edge.source.plot(index_mask = edge.matches['destination_idx'], alpha = 0)

        matplotlib.pyplot.scatter(kps['x'], kps['y'], color = 'black', alpha = 1)

```