Coupled Decomposition - Part III

import math
import os
import sys

sys.path.insert(0, os.path.abspath('/data/autocnet'))

import autocnet
from autocnet import CandidateGraph

# The GPU based extraction library that contains SIFT extraction and matching
import cudasift as cs

# A method to resize the images on the fly.
from scipy.misc import imresize

# Fundamental matrix computation
from autocnet.transformation import fundamental_matrix as fm

from autocnet.transformation.decompose import coupled_decomposition
from scipy.spatial.distance import cdist

%pylab inline
figsize(16,4)

Populating the interactive namespace from numpy and matplotlib

Preprocessing

a = 'AS15-P-0111_CENTER_LRG_CROPPED.png'
b = 'AS15-P-0112_CENTER_LRG_CROPPED.png'

adj = {a:[b],
       b:[a]}

cg = CandidateGraph.from_adjacency(adj)

# Enable the GPU
autocnet.cuda(enable=True, gpu=0)

# Write a custom keypoint extraction function - this could get monkey patched onto the graph object...
def extract(arr, downsample_amount=None, **kwargs):
    total_size = arr.shape[0] * arr.shape[1]
    if not downsample_amount:
        downsample_amount = math.ceil(total_size / 12500**2)
    shape = (int(arr.shape[0] / downsample_amount), int(arr.shape[1] / downsample_amount))
    # Downsample
    arr = imresize(arr, shape, interp='lanczos')

    npts = max(arr0.shape) / 3.5
    sd = cs.PySiftData(npts)
    cs.ExtractKeypoints(arr, sd, **kwargs)
    kp, des = sd.to_data_frame()
    kp = kp[['x', 'y', 'scale', 'sharpness', 'edgeness', 'orientation', 'score', 'ambiguity']]
    kp['score'] = 0.0
    kp['ambiguity'] = 0.0

    return kp, des, sd, downsample_amount, arr

# Write a generic decomposer
def custom_decompose(arr0, arr1):
    kp0, des0, sd0, downsample_amount0, arr0 = extract(arr0, thresh=1)
    kp1, des1, sd1, downsample_amount1, arr1 = extract(arr1, thresh=1)

    # Now apply matching, outlier detection, and compute a fundamental matrix
    sd0 = cs.PySiftData.from_data_frame(kp0, des0)
    sd1 = cs.PySiftData.from_data_frame(kp1, des1)

    # Apply the matcher
    cs.PyMatchSiftData(sd0, sd1)
    matches, _ = sd0.to_data_frame()
    # Generic decision about ambiguity and score based on quantiles

    # Apply outlier detection methods for the matches
    ambiguity_threshold = matches.ambiguity.quantile(0.01)  # Grabbing the 1%s in this data set
    score = matches.score.quantile(0.85)
    submatches = matches.query('ambiguity <= {} and score >= {}'.format(ambiguity_threshold, score))

    # Compute a fundamental matrix
    kpa = submatches[['x','y']]
    kpb = submatches[['match_xpos', 'match_ypos']]
    F, mask = fm.compute_fundamental_matrix(kpa, kpb, method='ransac', reproj_threshold=2.0)
    F = fm.enforce_singularity_constraint(F)

    # Grab the inliers
    inliers = submatches[mask]

    # Prepare for coupled decomposition
    midx = arr0.shape[1] / 2
    midy = arr0.shape[0] / 2

    mid = np.array([[midx, midy]])
    dists = cdist(mid, inliers[['x', 'y']])
    mid_correspondence = inliers.iloc[np.argmin(dists)]
    mid_correspondence

    # Decompose the images into quadrants
    smembership, dmembership, = coupled_decomposition(arr0, arr1,
                                                 sorigin=mid_correspondence[['x', 'y']],
                                                 dorigin=mid_correspondence[['match_xpos', 'match_ypos']],
                                                 theta=0)

    # Return the membership decisions
    return smembership, dmembership

import pandas as pd

gd0 = cg.node[0].geodata
gd1 = cg.node[1].geodata
membership0 = np.zeros(gd0.raster_size[::-1], dtype=np.int8)
membership1 = np.zeros(gd1.raster_size[::-1], dtype=np.int8)

# Recursively decompose twice.
pcounter = 0
for k in range(2):
    print(k)
    partitions = np.unique(membership0)
    for p in partitions:
        # Get the source extent as MBR
        sypart, sxpart = np.where(membership0 == p)
        minsy = np.min(sypart)
        maxsy = np.max(sypart) + 1
        minsx = np.min(sxpart)
        maxsx = np.max(sxpart) + 1

        del sypart, sxpart

        # Get the destination extent as MBR
        dypart, dxpart = np.where(membership1 == p)
        mindy = np.min(dypart)
        maxdy = np.max(dypart) + 1
        mindx = np.min(dxpart)
        maxdx = np.max(dxpart) + 1

        # Offsets into classified array
        sdy = dypart - min(dypart)
        sdx = dxpart - min(dxpart)

        print(minsy, maxsy, minsx, maxsx)
        print(mindy, maxdy, mindx, maxdx)

        arr0 = gd0.read_array(pixels=map(int, [minsx, minsy, maxsx-minsx, maxsy-minsy]))
        arr1 = gd1.read_array(pixels=map(int,[mindx, mindy, maxdx-mindx, maxdy-mindy]))

        smem, dmem = custom_decompose(arr0, arr1)
        smem += pcounter
        dmem += pcounter

        # Some fancy indexing to get the dmembership into the right place in the global membership
        membership0[minsy:maxsy,minsx:maxsx] = smem
        membership1[dypart, dxpart] = dmem[sdy, sdx]
        pcounter += 4

        # Force cleanup
        arr0 = None
        arr1 = None

0
0 11510 0 59720
0 11568 0 59770
0 0 59720 11510
0 0 59770 11568

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-a10dc6618433> in <module>()
     43
     44         # Some fancy indexing to get the dmembership into the right place in the global membership
---> 45         membership0[minsy:maxsy,minsx:maxsx] = smem
     46         membership1[dypart, dxpart] = dmem[sdy, sdx]
     47         pcounter += 4

ValueError: could not broadcast input array from shape (2302,11944) into shape (11510,59720)