Added more voronoi notebook code to the code base.

import pandas as pd

import numpy as np

import networkx as nx

import geopandas as gpd

import ogr

    return intersection_area

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

def vor(graph, clean_keys, s=30):

        """

        Creates a voronoi diagram for an edge using either the coordinate

        transformation or using the homography between source and destination.

        Creates a voronoi diagram for a complete graph using either the coordinate

        transformation or using the homography based around a source node.

        The coordinate transformation uses the footprint of source and destination to

        calculate an intersection between the two images, then transforms the vertices of

        The coordinate transformation uses the footprint of each image to

        calculate an intersection between the image set, then transforms the vertices of

        the intersection back into pixel space.

        If a coordinate transform does not exist, use the homography to project the destination image

        If a coordinate transform does not exist, use the homography to project the n - 1 images

        onto the source image, producing an area of intersection.

        The intersection vertices are then scaled by a factor of s (default 30), this accounts for the

        Parameters

        ----------

        edge : object

               An edge object

        graph : object

               A networkx graph object

        clean_keys : list

                     Of strings used to apply masks to omit correspondences

                     3 column pandas dataframe of x, y, and weights

        """

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

        for n in graph.nodes():

            neighbors = len(graph.neighbors(n))

            if neighbors != num_neighbors:

        neighbors_dict = nx.degree(graph)

        if not all(value == len(graph.neighbors(graph.nodes()[0])) for value in neighbors_dict.values()):

                raise AssertionError('The graph is not complete')

        intersection, proj_gdf, source_gdf = compute_intersection(graph.nodes()[0], 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)

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

                                       'weight'] = poly_area

                i += 1

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

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

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

    source_corners = source.geodata.xy_corners

    source_poly = Polygon(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

        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:

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

        if (source.geodata.coordinate_transformation.this is None) and \

           (destination.geodata.coordinate_transformation.this is None):

            proj_poly = Polygon(destination.geodata.xy_corners)

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

from plio.io import io_hdf, io_json

from plio.utils import utils as io_utils

from plio.io.io_gdal import GeoDataset

from autocnet.cg.cg import vor

from autocnet.graph import markov_cluster

from autocnet.graph.edge import Edge

from autocnet.graph.node import Node

        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

    def compute_cliques(self, node_id):

        if node_id is not None:

            return list(nx.cliques_containing_node(self, nodes=node_id))

        else:

            return list(nx.find_cliques(self))

    def compute_vor_weight(self, clean_keys, node_id=None, clique=False):

        if clique:

            cliques = self.compute_cliques(node_id)

            for g in cliques:

                subgraph = self.create_node_subgraph(g)

                vor(subgraph, clean_keys)

        else:

            for e in self.edges_iter():

                subgraph = self.create_node_subgraph(e)

                vor(subgraph, clean_keys)

import numpy as np

import pandas as pd

import networkx as nx

from osgeo import ogr

    centroid = centroid[:2]

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

    return vector

%% 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.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')

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

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

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

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

    neighbors_dict = nx.degree(graph)

    if not all(value == len(graph.neighbors(graph.nodes()[0])) for value in neighbors_dict.values()):

        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['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()

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

def compute_weights(cang, clean_keys = [], node_id=None, clique = False,):

    if clique:

        cliques = compute_cliques(cang, node_id)

        for g in cliques:

            subgraph = cang.create_node_subgraph(g)

            cg.vor(subgraph, clean_keys)

    else:

        for e in cang.edges_iter():

            subgraph = cang.create_node_subgraph(e)

            cg.vor(subgraph, clean_keys)

```

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

def compute_cliques(graph, node_id):

    if node_id is not None:

        return list(nx.cliques_containing_node(graph, nodes=node_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'])

compute_weights(cang, clique = True, clean_keys = ['fundamental'])

print(time.clock() - start)

start = time.clock()

compute_weights(cang, clique = False, clean_keys = ['fundamental'])

print(time.clock() - start)

```

%% Output

    3.1514979999999966

    2.8140520000000038

%% Cell type:code id: tags:

``` python

start = time.clock()

cang.compute_vor_weight(clean_keys = ['fundamental'], clique = True)

print(time.clock() - start)

start = time.clock()

for e in cang.edges_iter():

    subgraph = cang.create_node_subgraph(e)

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

cang.compute_vor_weight(clean_keys = ['fundamental'], clique = False)

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 >

    2.850580000000001

    2.9477230000000016

    3.1148349999999994

    2.8084049999999934

%% 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.__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