Subpixel (Part I - Changes to the underlying API) (#364)

    apriorisample = Column(Float)

    samplesigma = Column(Float)

    linesigma = Column(Float)

    weight = Column(Float, default=1)

    rms = Column(Float)

    @hybrid_property

from scipy.ndimage.interpolation import zoom

def pattern_match(template, image, upsampling=16, func=cv2.TM_CCOEFF_NORMED, error_check=False):

def pattern_match(template, image, upsampling=16, func=cv2.TM_CCORR_NORMED, error_check=False):

    """

    Call an arbitrary pattern matcher

               The strength of the correlation in the range [-1, 1].

    """

    different = {cv2.TM_SQDIFF_NORMED: cv2.TM_CCOEFF_NORMED,

                 cv2.TM_CCORR_NORMED: cv2.TM_SQDIFF_NORMED,

                 cv2.TM_CCOEFF_NORMED: cv2.TM_SQDIFF_NORMED}

    if upsampling < 1:

        raise ValueError

    # Fit a 3rd order polynomial to upsample the images

    u_template = zoom(template, upsampling, order=3)

    u_image = zoom(image, upsampling, order=3)

    result = cv2.matchTemplate(u_image, u_template, method=func)

    min_corr, max_corr, min_loc, max_loc = cv2.minMaxLoc(result)

    _, max_corr, min_loc, max_loc = cv2.minMaxLoc(result)

    if func == cv2.TM_SQDIFF or func == cv2.TM_SQDIFF_NORMED:

        x, y = (min_loc[0], min_loc[1])

    else:

import json

from math import modf, floor

import numpy as np

from skimage.feature import register_translation

from redis import StrictRedis

from plurmy import Slurm

from autocnet import Session, config

from autocnet.matcher import naive_template

from autocnet.matcher import ciratefi

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

from autocnet.graph.node import NetworkNode

import geopandas as gpd

import pandas as pd

            (nmatches, 1) to store the updated y coordinates

    """

    # Setup to store output to append to dataframes

    shifts_x = np.empty(nmatches)

    shifts_x[:] = np.nan

    shifts_y = np.empty(nmatches)

    shifts_y[:] = np.nan

    strengths = np.empty((nmatches, nstrengths))

    strengths[:] = np.nan

    shifts_x = np.zeros(nmatches)

    shifts_y = np.zeros(nmatches)

    strengths = np.zeros((nmatches, nstrengths))

    new_x = np.empty(nmatches)

    new_y = np.empty(nmatches)

    return shifts_x, shifts_y, strengths, new_x, new_y

def check_image_size(imagesize):

    """

    Given an x,y tuple, ensure that the values

    are odd. Used by the subpixel template to also ensure

    that the template size is the one requested and not 2x

    the template size.

    Parameters

    ----------

    imagesize : tuple

                in the form (size_x, size_y)

    """

    x = imagesize[0] / 2

    y = imagesize[1] / 2

    if x % 2 == 0:

        x += 1

    if y % 2 == 0:

        y += 1

    return x,y

def clip_roi(img, center_x, center_y, size_x=200, size_y=200):

    """

    Given an input image, clip a square region of interest

          with a read_array method that takes a pixels

          argument in the form [xstart, ystart, xstop, ystop]

    center : tuple

             (x,y) coordinates to center the roi

    center_x : Numeric

               The x coordinate to the center of the roi

    center_y : Numeric

               The y coordinate to the center of the roi

    img_size : int

               Odd, total image size

               1/2 of the total image size. This value is the

               number of pixels grabbed from each side of the center

    Returns

    -------

        size_y = int(ay)

    # Read from the upper left origin

    pixels=(int(ax-size_x), int(ay-size_y), size_x * 2, size_y * 2)

    pixels = [ax-size_x, ay-size_y, size_x*2, size_y*2]

    pixels = list(map(int, pixels))  # 

    if isinstance(img, np.ndarray):

        subarray = img[pixels[1]:pixels[1] + pixels[3] + 1, pixels[0]:pixels[0] + pixels[2] + 1]

    else:

    search : ndarray

             The search image

    Returns

    Returnsslurm-2235260_89.out.2235260_89.out

    -------

    x_offset : float

               Shift in the x-dimension

    (y_shift, x_shift), error, diffphase = register_translation(search, template, **kwargs)

    return x_shift, y_shift, (error, diffphase)

def subpixel_template(sx, sy, dx, dy, s_img, d_img, search_size=251, template_size=51, **kwargs):

def subpixel_template(sx, sy, dx, dy, s_img, d_img, image_size=(251, 215), template_size=(51,51), **kwargs):

    """

    Uses a pattern-matcher on subsets of two images determined from the passed-in keypoints and optional sizes to

    compute an x and y offset from the search keypoint to the template keypoint and an associated strength.

    Parameters

    ----------

    sx : Numeric

         Source X coordinate

    sy : Numeric

         Source y coordinate

    dx : Numeric

         The desintation x coordinate

    dy : Numeric

         The destination y coordinate

    s_img : GeoDataset

            The source image GeoDataset

    d_img : GeoDataset

            The destination image GeoDataset

    image_size : tuple

                 (xsize, ysize) of the image that is searched within (this should be larger 

                 than the template size)

    template_size : tuple

                    (xsize, ysize) of the template to iterate over the image in order

                    to identify the area(s) of highest correlation.

    Returns

    -------

    x_shift : float

               Shift in the x-dimension

    y_shift : float

               Shift in the y-dimension

    strength : float

               Strength of the correspondence in the range [-1, 1]

    See Also

    --------

    autocnet.matcher.naive_template.pattern_match : for the kwargs that can be passed to the matcher

    """

    image_size = check_image_size(image_size)

    template_size = check_image_size(template_size)

    s_image, _, _ = clip_roi(s_img, sx, sy, size_x=image_size[0], size_y=image_size[1])

    d_template, dxr, dyr = clip_roi(d_img, dx, dy, size_x=template_size[0], size_y=template_size[1])

    if (s_image is None) or (d_template is None):

        return None, None, None

    shift_x, shift_y, metrics = naive_template.pattern_match(d_template, s_image, **kwargs)

    dx = (dx - shift_x + dxr)

    dy = (dy - shift_y + dyr)  

    return dx, dy, metrics

    x_shift, y_shift, strength = naive_template.pattern_match(template, image, **kwargs)

    return x_shift, y_shift, strength

def subpixel_ciratefi(sx, sy, dx, dy, s_img, d_img, search_size=251, template_size=51, **kwargs):

    """

    Uses a pattern-matcher on subsets of two images determined from the passed-in keypoints and optional sizes to

    compute an x and y offset from the search keypoint to the template keypoint and an associated strength.

    strength : float

               Strength of the correspondence in the range [-1, 1]

    """

    template, _, _ = clip_roi(d_img, dx, dy,

                              size_x=template_size, size_y=template_size)

    search, dxr, dyr = clip_roi(s_img, sx, sy,

                                size_x=search_size, size_y=search_size)

    if template is None or search is None:

        return None, None, None

    if 'method' in kwargs.keys():

        method = kwargs['method']

        kwargs.pop('method', None)

    else:

        method = 'naive'

    functions = { 'naive' : naive_template.pattern_match,

                  'ciratefi' : ciratefi.ciratefi}

    x_offset, y_offset, strength = functions[method](template, search, **kwargs)

    x_offset, y_offset, strength = ciratefi.ciratefi(template, search, **kwargs)

    dx += (x_offset + dxr)

    dy += (y_offset + dyr)

    return dx, dy, strength

def iterative_phase(sx, sy, dx, dy, s_img, d_img, size=251, reduction=11, convergence_threshold=1.0, max_dist=50, **kwargs):

    """

    Iteratively apply a subpixel phase matcher to source (s_img) amd destination (d_img)

    Iteratively apply a subpixel phase matcher to source (s_img) and destination (d_img)

    images. The size parameter is used to set the initial search space. The algorithm

    is recursively applied to reduce the total search space by reduction until the convergence criteria

    are met. Convergence is defined as the point at which the computed shifts (x_shift,y_shift) are

            break

    return dx, dy, metrics

def subpixel_register_point(pointid, iterative_phase_kwargs={}, subpixel_template_kwargs={},

                            cost_func=lambda x,y: 1/x**2 * y, threshold=0.005):

    """

    Given some point, subpixel register all of the measures in the point to the 

    first measure.

    Parameters

    ----------

    pointid : int

              The identifier of the point in the DB

    iterative_phase_kwargs : dict

                             Any keyword arguments passed to the phase matcher

    subpixel_template_kwargs : dict

                               Ay keyword arguments passed to the template matcher

    cost : func

           A generic cost function accepting two arguments (x,y), where x is the 

           distance that a point has shifted from the original, sensor identified 

           intersection, and y is the correlation coefficient coming out of the 

           template matcher.

    threshold : numeric

                measures with a cost <= the threshold are marked as active=False in

                the database.

    """

    session = Session()

    point = session.query(Points).filter(Points.id == pointid).one()

    measures = point.measures

    source = measures[0]

    sourceid = source.imageid

    res = session.query(Images).filter(Images.id == sourceid).one()

    source_node = NetworkNode(node_id=sourceid, image_path=res.path)

    for measure in measures[1:]:

        active = True

        cost = None

        destinationid = measure.imageid

        res = session.query(Images).filter(Images.id == destinationid).one()

        destination_node = NetworkNode(node_id=destinationid, image_path=res.path)

        new_phase_x, new_phase_y, phase_metrics = iterative_phase(source.sample, 

                                                                source.line, 

                                                                measure.sample, 

                                                                measure.line,

                                                                source_node.geodata, 

                                                                destination_node.geodata,

                                                                **iterative_phase_kwargs)

        if new_phase_x == None:

            measure.active = False # Unable to phase match

            continue

        new_template_x, new_template_y, template_metric = subpixel_template(source.sample,

                                                                source.line,

                                                                new_phase_x,

                                                                new_phase_y,

                                                                source_node.geodata,

                                                                destination_node.geodata,

                                                                **subpixel_template_kwargs)

        if new_template_x == None:

            measure.active = False # Unable to template match

            continue

        dist = np.linalg.norm([new_phase_x-new_template_x, new_phase_y-new_template_y])

        cost = cost_func(dist, template_metric)

        if cost <= threshold:

            measure.active = False # Threshold criteria not met

            continue

        # Update the measure

        measure.active = active

        if new_template_x:

            measure.sample = new_template_x

            measure.line = new_template_y

            measure.weight = cost

    session.commit()

    session.close()

def subpixel_register_points(iterative_phase_kwargs={'size': 71}, 

                             subpixel_template_kwargs={'image_size':(121,121)},

                             cost_func=lambda x,y: 1/x**2 * y,

                             threshold=0.005):

    """

    Serial subpixel registration of all of the points in a given DB table.

    Parameters

    ----------

    pointid : int

              The identifier of the point in the DB

    iterative_phase_kwargs : dict

                             Any keyword arguments passed to the phase matcher

    subpixel_template_kwargs : dict

                               Ay keyword arguments passed to the template matcher

    cost : func

           A generic cost function accepting two arguments (x,y), where x is the 

           distance that a point has shifted from the original, sensor identified 

           intersection, and y is the correlation coefficient coming out of the 

           template matcher.

    threshold : numeric

                measures with a cost <= the threshold are marked as active=False in

                the database.

    """

    session = Session()

    pointids = [point.id for point in session.query(Points)]

    sessoin.close()

    for pointid in pointids:

        subpixel_register_point(pointid, 

                                iterative_phase_kwargs=iterative_phase_kwargs,

                                subpixel_template_kwargs=subpixel_template_kwargs, 

                                cost_func=cost_func)

def cluster_subpixel_register_points(iterative_phase_kwargs={'size': 71}, 

                                     subpixel_template_kwargs={'image_size':(121,121)},

                                     cost_func=lambda x,y: 1/x**2 * y,

                                     threshold=0.005,

                                     walltime='00:10:00'):

    """

    Distributed subpixel registration of all of the points in a given DB table.

    Parameters

    ----------

    pointid : int

              The identifier of the point in the DB

    iterative_phase_kwargs : dict

                             Any keyword arguments passed to the phase matcher

    subpixel_template_kwargs : dict

                               Ay keyword arguments passed to the template matcher

    cost : func

           A generic cost function accepting two arguments (x,y), where x is the 

           distance that a point has shifted from the original, sensor identified 

           intersection, and y is the correlation coefficient coming out of the 

           template matcher.

    threshold : numeric

                measures with a cost <= the threshold are marked as active=False in

                the database.

    """    

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

    session = Session()

    for i, point in enumerate(session.query(Points)):

        msg = {'id' : point.id,

               'iterative_phase_kwargs' : iterative_phase_kwargs,

               'subpixel_template_kwargs' : subpixel_template_kwargs,

               'threshold':threshold,

               'walltime' : walltime}

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

    session.close()

    job_counter = i + 1

    # Submit the jobs

    submitter = Slurm('acn_subpixel',

                 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

    arrs = sp._prep_subpixel(nmatches, nstrengths=nstrengths)

    assert len(arrs) == 5

    assert arrs[2].shape == (nmatches, nstrengths)

    assert np.isnan(arrs[0][0])

    assert arrs[0][0] == 0

@pytest.mark.parametrize("center_x, center_y, size, expected", [(4, 4, 9, 404),

                                                          (55.4, 63.1, 27, 6355)])

                                                a, b, upsampling=16)

    assert strength >= 0.99

    assert nx == 50.5625 

    assert ny == 52.5625 

    assert nx == 50.9375

    assert ny == 48.9375

@pytest.mark.parametrize("convergence_threshold, expected", [(1.0, (None, None, None)),

                                                             (2.0, (50.49, 52.44, (0.039507, -9.5e-20)))])

        if expected[2] is not None:

            for i in range(len(strength)):

                assert pytest.approx(strength[i],6) == expected[2][i]

@pytest.mark.parametrize("data, expected", [

    ((21,21), (10.5, 10.5)),

    ((20,20), (11,11)),

    ((0,0), (1,1))

])

def test_check_image_size(data, expected):

    assert sp.check_image_size(data) == expected

 No newline at end of file

#!/usr/bin/env python

import copy

import os

import json

import sys

from redis import StrictRedis

import yaml

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

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

from autocnet.matcher.subpixel import subpixel_register_point

from autocnet.graph.node import NetworkNode

from autocnet import Session, dem

from plio.io.io_gdal import GeoDataset

#Load the config file

try:

    with open(os.environ['autocnet_config'], 'r') as f:

        config = yaml.safe_load(f)

except:

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

    sys.exit(1)

def main(msg, config):

    session = Session()

    id = msg['id']

    point = session.query(Points).filter(Points.id == id).one()

    pointid = point.id

    session.close()

    if point is None:

        print(f'Unable to find point {pointid}')

    print(f'Subpixel registering point {pointid}')

    subpixel_register_point(pointid,

                            iterative_phase_kwargs=msg['iterative_phase_kwargs'],

                            subpixel_template_kwargs=msg['subpixel_template_kwargs'],

                            threshold=msg['threshold'])

    print('Subpixel registration complete.')

if __name__ == '__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)

 No newline at end of file