Commit 387cc47a authored by Adam Paquette's avatar Adam Paquette
Browse files

Seperated some voronoi code into helper methods that may have uses elsewhere....

Seperated some voronoi code into helper methods that may have uses elsewhere. Also started adding notebook code to the code base.
parent 728affd7
Loading
Loading
Loading
Loading
+48 −38
Changes for autocnet/cg/cg.py: 48 added lines, 38 removed lines.
Original line number Diff line number Diff line
@@ -8,9 +8,10 @@ import ogr

from scipy.spatial import ConvexHull
from scipy.spatial import Voronoi
import shapely.geometry
from shapely.geometry import Polygon, Point
from shapely.affinity import scale
from shapely.ops import unary_union
# from shapely.ops import unary_union
import cv2

from autocnet.utils import utils
@@ -61,6 +62,19 @@ def convex_hull(points):
    return hull


def geom_mask(keypoints, geom):
    def _in_mbr(r, mbr):
        if (mbr[0] <= r.x <= mbr[2]) and (mbr[1] <= r.y <= mbr[3]):
            return True
        else:
            return False

    mbr = geom.bounds
    initial_mask = keypoints.apply(_in_mbr, axis=1, args=(mbr,))

    return initial_mask


def two_poly_overlap(poly1, poly2):
    """

@@ -114,7 +128,7 @@ def get_area(poly1, poly2):
    return intersection_area


def vor(graph, clean_keys, s=30):
def compute_voronoi(keypoints, intersection=None, geometry=False, s=30):
        """
        Creates a voronoi diagram for all edges in a graph, and assigns a given
        weight to each edge. This is based around voronoi polygons generated
@@ -131,59 +145,55 @@ def vor(graph, clean_keys, s=30):
        s : int
            Offset for the corners of the image
        """
        neighbors_dict = nx.degree(graph)
        if not all(value == len(graph.neighbors(graph.nodes()[0])) for value in neighbors_dict.values()):
            warnings.warn('The given graph is not complete and may yield garbage.')

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

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

        source_node = graph.nodes()[0]
        for e in graph.edges():
            # If the source node changes, change the projection space to match it
            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]]
            kps = edge.get_keypoints('source', clean_keys=clean_keys, homogeneous=True)
        vor_keypoints = []

            kps['geometry'] = kps.apply(lambda x: Point(x['x'], x['y']), axis=1)
            kps.mask = kps['geometry'].apply(lambda x: intersection.geometry.contains(x).any())
        keypoints.apply(lambda x: vor_keypoints.append((x['x'], x['y'])), axis = 1)

            # Creates a mask for displaying the voronoi points
            # Currently erronious and produces NaN values in place
            # of true
            # matches, mask = edge.clean(clean_keys = clean_keys)
            # mask[mask] = kps.mask
            # edge.masks = ('voronoi', mask)
        if intersection is None:
            keypoint_bounds = Polygon(vor_keypoints).bounds
            min_bounding_box = shapely.geometry.box(keypoint_bounds[0], keypoint_bounds[1],
                                                    keypoint_bounds[2], keypoint_bounds[3])

            keypoints = []
            kps[kps.mask].apply(lambda x: keypoints.append((x['x'], x['y'])), axis=1)
            scaled_coords = np.array(scale(min_bounding_box, s, s).exterior.coords)
        else:
            scaled_coords = np.array(scale(intersection, s, s).exterior.coords)

            scaled_coords = np.array(scale(source_gdf.geometry[0], s, s).exterior.coords)
            keypoints = np.vstack((keypoints, scaled_coords))
        vor_keypoints = np.vstack((vor_keypoints, scaled_coords))
        vor = Voronoi(vor_keypoints)

            vor = Voronoi(keypoints)
            voronoi_df = pd.DataFrame(data=kps, columns=['x', 'y', 'weight'])
        # For weight computation
        # Should move to its own method
        if geometry:
            voronoi_df = gpd.GeoDataFrame(data = keypoints, columns=['x', 'y', 'weight', 'geometry'])
        else:
            voronoi_df = gpd.GeoDataFrame(data = keypoints, 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 -1 not in region:

            if not -1 in region:
                polygon_points = [vor.vertices[i] for i in region]

                if len(polygon_points) != 0:
                    polygon = Polygon(polygon_points)
                        poly_area = polygon.intersection(intersection.geometry[intersection_ind]).area

                    if intersection is not None:
                        intersection_poly = polygon.intersection(intersection)
                    else:
                        intersection_poly = polygon.intersection(min_bounding_box)

                    voronoi_df.loc[(voronoi_df["x"] == region_point[0][0][0]) &
                                   (voronoi_df["y"] == region_point[0][0][1]),
                                   'weight'] = intersection_poly.area
                    if geometry:
                        voronoi_df.loc[(voronoi_df["x"] == region_point[0][0][0]) &
                                       (voronoi_df["y"] == region_point[0][0][1]),
                                       'weight'] = poly_area
                                       'geometry'] = intersection_poly
            i += 1

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


def compute_intersection(graph, source, clean_keys=[]):
+4 −1
Changes for autocnet/graph/edge.py: 4 added lines, 1 removed line.
Original line number Diff line number Diff line
@@ -454,7 +454,7 @@ class Edge(dict, MutableMapping):
        """
        pass

    def get_keypoints(self, node, clean_keys):
    def get_keypoints(self, node, clean_keys, homogeneous=False):
        """

        Returns a list of keypoint coordinates that match the specified
@@ -514,4 +514,7 @@ class Edge(dict, MutableMapping):
        # Return keypts @ masked indecies for the node
        masked_keypts = all_keypts.iloc[keypt_indices].sort_index()

        if homogeneous:
            masked_keypts['homogeneous'] = 1

        return masked_keypts
+79 −3
Changes for autocnet/graph/network.py: 79 added lines, 3 removed lines.
Original line number Diff line number Diff line
@@ -5,11 +5,13 @@ import warnings

import networkx as nx
import pandas as pd
import shapely.geometry

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.cg.cg import geom_mask
from autocnet.cg.cg import compute_voronoi
from autocnet.graph import markov_cluster
from autocnet.graph.edge import Edge
from autocnet.graph.node import Node
@@ -706,7 +708,7 @@ class CandidateGraph(nx.Graph):
        else:
            return list(nx.find_cliques(self))

    def compute_vor_weight(self, clean_keys):
    def compute_weight(self, clean_keys):
        """
        Computes a voronoi weight for each edge in a given graph.
        Can function as is, but is slightly optimized for complete subgraphs.
@@ -715,6 +717,80 @@ class CandidateGraph(nx.Graph):
        ----------
        clean_keys : list
                     Strings used to apply masks to omit correspondences
        """
        neighbors_dict = nx.degree(self)
        if False in list(all(value == len(self.neighbors(self.nodes()[0])) for value in neighbors_dict.values())):
            warnings.warn('The given graph is not complete and may yield garbage.')

        intersect_gdf = self.compute_intersection(self.nodes()[0], self, clean_keys)
        source_node = self.node[self.nodes()[0]]

        for s, d, edge in self.edges_iter(data=True):
            # Recompute the intersection if the source node of the n + 1 edge is different from the n edge
            # I don't know if this check is necessary anymore
            if s != source_node['node_id']:
                source_node = edge.source
                intersect_gdf = self.compute_intersection(self, source_node, clean_keys)

            kps = edge.get_keypoints('source', clean_keys=clean_keys)[['x', 'y']]
            reproj_geom = source_node.reproject_geom(intersect_gdf.query("overlaps_all == True").geometry.values[0].__geo_interface__['coordinates'][0])
            initial_mask = geom_mask(kps, reproj_geom)

            if (len(kps[initial_mask]) <= 0):
                continue

            kps['geometry'] = kps.apply(lambda x: shapely.geometry.Point(x['x'], x['y']), axis=1)
            kps_mask = kps['geometry'][initial_mask].apply(lambda x: reproj_geom.contains(x))
            voronoi_df = compute_voronoi(kps[initial_mask][kps_mask], reproj_geom, geometry=True)

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

    def compute_intersection(self, source, clean_keys=[], ax = None):
        """
        Computes the intercetion of all images in a graph
        based around a given source node

        Parameters
        ----------
        source: object or int
                    Either a networkx Node object or an integer

        clean_keys : list
                     Strings used to apply masks to omit correspondences
        """
        vor(self, clean_keys)
        if type(source) is int:
            source = self.node[source]

        try:
            source_poly = swkt.loads(source.geodata.footprint.GetGeometryRef(0).ExportToWkt())
        except:
            raise AttributeError()

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

        proj_list = []
        proj_nodes = []

        # Begin iterating through the nodes in the graph excluding the source node
        for s, d, edge in self.edges_iter(data=True):
            if s == source['node_id']:
                proj_poly = swkt.loads(edge.destination.geodata.footprint.GetGeometryRef(0).ExportToWkt())
            elif d == source['node_id']:
                proj_poly = swkt.loads(edge.source.geodata.footprint.GetGeometryRef(0).ExportToWkt())
            else:
                continue

            proj_list.append(proj_poly)
            proj_nodes.append(n)


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

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

        if len(intersect_gdf.query("overlaps_all == True")) <= 0:
            new_poly = unary_union(intersect_gdf.geometry)
            intersect_gdf.loc[len(intersect_gdf)] = [source['node_id'], source['node_id'], new_poly, True]

        return intersect_gdf
+8 −0
Changes for autocnet/graph/node.py: 8 added lines, 0 removed lines.
Original line number Diff line number Diff line
@@ -7,6 +7,7 @@ import pandas as pd
from plio.io.io_gdal import GeoDataset
from plio.io.isis_serial_number import generate_serial_number
from scipy.misc import bytescale
import shapley

from autocnet.cg import cg
from autocnet.control.control import Correspondence, Point
@@ -511,3 +512,10 @@ class Node(dict, MutableMapping):
        mask = panel[clean_keys].all(axis=1)
        matches = self._keypoints[mask]
        return matches, mask

    def reproject_geom(self, coords):
        reproj = []

        for x, y in coords:
            reproj.append(self.geodata.latlon_to_pixel(x, y))
        return shapley.Polygon(reproj)
+470 −409

File changed.

Preview size limit exceeded, changes collapsed.