Commit 6fd6b2b6 authored by Paquette, Adam Chance's avatar Paquette, Adam Chance
Browse files

Merge branch 'api' into 'main'

Removing ciratefi, add MI matcher from subpixel api

See merge request astrogeology/autocnet!663
parents 27705fe7 3c4e5ba5
Loading
Loading
Loading
Loading

autocnet/matcher/ciratefi.py

deleted100644 → 0
+0 −758

File deleted.

Preview size limit exceeded, changes collapsed.

+17 −103
Original line number Diff line number Diff line
from math import floor
import logging

from autocnet.transformation.roi import Roi
import numpy as np

from scipy.ndimage.measurements import center_of_mass
import skimage.transform as tf

import logging 
# setup logging file
log = logging.getLogger(__name__)

def mutual_information(t1, t2, **kwargs):

def mutual_information(reference_arr, moving_arr, **kwargs):
    """
    Computes the correlation coefficient between two images using a histogram
    comparison (Mutual information for joint histograms). The corr_map coefficient
@@ -17,12 +19,13 @@ def mutual_information(t1, t2, **kwargs):
    Parameters
    ----------

    t1 : ndarray
    reference_arr : ndarray
                    First image to use in the histogram comparison
    
    t2 : ndarray
    moving_arr : ndarray
                   Second image to use in the histogram comparison
    
    
    Returns
    -------

@@ -35,15 +38,15 @@ def mutual_information(t1, t2, **kwargs):
    numpy.histogram2d : for the kwargs that can be passed to the comparison
    """
    
    if np.isnan(t1).any() or np.isnan(t2).any():
        log.warning('Unable to process due to NaN values in the input data')
    if np.isnan(reference_arr).any() or np.isnan(moving_arr).any():
        log.warning('Unable compute MI. Image sizes are not identical.')        
        return
    
    if t1.shape != t2.shape:
    if reference_arr.shape != moving_arr.shape:
        log.warning('Unable compute MI. Image sizes are not identical.')
        return

    hgram, x_edges, y_edges = np.histogram2d(t1.ravel(),t2.ravel(), **kwargs)
    hgram, x_edges, y_edges = np.histogram2d(reference_arr.ravel(),moving_arr.ravel(), **kwargs)

    # Convert bins counts to probability values
    pxy = hgram / float(np.sum(hgram))
@@ -54,92 +57,3 @@ def mutual_information(t1, t2, **kwargs):
    nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
    return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
def mutual_information_match(d_template, s_image, subpixel_size=3,
                             func=None, **kwargs):
    """
    Applys the mutual information matcher function over a search image using a
    defined template


    Parameters
    ----------
    d_template : ndarray
                 The input search template used to 'query' the destination
                 image

    s_image : ndarray
              The image or sub-image to be searched

    subpixel_size : int
                    Subpixel area size to search for the center of mass
                    calculation

    func : function
           Function object to be used to compute the histogram comparison

    Returns
    -------
    x : float
        The x offset

    y : float
        The y offset

    max_corr : float
               The strength of the correlation in the range [0, 4].

    corr_map : ndarray
               Map of corrilation coefficients when comparing the template to
               locations within the search area
    """

    if func == None:
        func = mutual_information

    image_size = s_image.shape
    template_size = d_template.shape

    y_diff = image_size[0] - template_size[0]
    x_diff = image_size[1] - template_size[1]

    max_corr = -np.inf
    corr_map = np.zeros((y_diff+1, x_diff+1))
    max_i = -1  # y
    max_j = -1  # x
    for i in range(y_diff+1):
        for j in range(x_diff+1):
            sub_image = s_image[i:i+template_size[1],  # y
                                j:j+template_size[0]]  # x
            corr = func(sub_image, d_template, **kwargs)
            if corr > max_corr:
                max_corr = corr
                max_i = i
                max_j = j
            corr_map[i, j] = corr

    y, x = np.unravel_index(np.argmax(corr_map, axis=None), corr_map.shape)

    upper = int(2 + floor(subpixel_size / 2))
    lower = upper - 1
    # x, y are the location of the upper left hand corner of the template in the image
    area = corr_map[y-lower:y+upper,
                    x-lower:x+upper]

    # Compute the y, x shift (subpixel) using scipys center_of_mass function
    cmass  = center_of_mass(area)

    if area.shape != (subpixel_size+2, subpixel_size+2):
        log.warning("Max correlation is too close to the boundary.")
        return None, None, 0, None

    subpixel_y_shift = subpixel_size - 1 - cmass[0]
    subpixel_x_shift = subpixel_size - 1 - cmass[1]

    y += subpixel_y_shift
    x += subpixel_x_shift

    # Compute the idealized shift (image center)
    y -= (s_image.shape[0] / 2) - (d_template.shape[0] / 2)
    x -= (s_image.shape[1] / 2) - (d_template.shape[1] / 2)

    return float(x), float(y), float(max_corr), corr_map
+89 −1
Original line number Diff line number Diff line
@@ -10,6 +10,7 @@ import logging

import numpy as np

from scipy.ndimage.measurements import center_of_mass
from skimage import transform as tf
from skimage import registration
from skimage import filters
@@ -28,7 +29,7 @@ import PIL
from PIL import Image

from autocnet.matcher.naive_template import pattern_match
from autocnet.matcher.mutual_information import mutual_information_match
from autocnet.matcher.mutual_information import mutual_information
from autocnet.spatial import isis
from autocnet.io.db.model import Measures, Points, Images, JsonEncoder
from autocnet.graph.node import NetworkNode
@@ -1400,3 +1401,90 @@ def smart_register_point(pointid, parameters=[], shared_kwargs={}, valid_reproje
    log.info(f'Ignoring measures: {measures_to_set_false}')

    return measures_to_update, measures_to_set_false


def mutual_information_match(moving_roi,
                             reference_roi, 
                             affine=tf.AffineTransform(), 
                             subpixel_size=3,
                             func=None, **kwargs):
    """
    Applies the mutual information matcher function over a search image using a
    defined template


    Parameters
    ----------
    moving_roi : roi 
                 The input search template used to 'query' the destination
                 image

    reference_roi : roi
              The image or sub-image to be searched

    subpixel_size : int
                    Subpixel area size to search for the center of mass
                    calculation

    func : function
           Function object to be used to compute the histogram comparison

    Returns
    -------
    new_affine :AffineTransform
                The affine transformation

    max_corr : float
               The strength of the correlation in the range [0, 4].

    corr_map : ndarray
               Map of corrilation coefficients when comparing the template to
               locations within the search area
    """
    reference_roi.clip()
    moving_roi.clip(affine=affine)

    moving_image = moving_roi.clipped_array
    reference_template = reference_roi.clipped_array

    if func == None:
        func = mutual_information

    image_size = moving_image.shape
    template_size = reference_template.shape

    y_diff = image_size[0] - template_size[0]
    x_diff = image_size[1] - template_size[1]

    max_corr = -np.inf
    corr_map = np.zeros((y_diff+1, x_diff+1))
    for i in range(y_diff+1):
        for j in range(x_diff+1):
            sub_image = moving_image[i:i+template_size[1],  # y
                                j:j+template_size[0]]  # x
            corr = func(sub_image, reference_template, **kwargs)
            if corr > max_corr:
                max_corr = corr
            corr_map[i, j] = corr

    y, x = np.unravel_index(np.argmax(corr_map, axis=None), corr_map.shape)

    upper = int(2 + floor(subpixel_size / 2))
    lower = upper - 1
    area = corr_map[y-lower:y+upper,
                    x-lower:x+upper]

    # Compute the y, x shift (subpixel) using scipys center_of_mass function
    cmass  = center_of_mass(area)
    if area.shape != (subpixel_size + 2, subpixel_size + 2):
        return  None, 0, None
        

    subpixel_y_shift = subpixel_size - 1 - cmass[0]
    subpixel_x_shift = subpixel_size - 1 - cmass[1]
    y = abs(y - (corr_map.shape[1])/2)
    x = abs(x - (corr_map.shape[0])/2)
    y += subpixel_y_shift
    x += subpixel_x_shift
    new_affine = tf.AffineTransform(translation=(-x, -y))
    return new_affine, np.max(max_corr), corr_map
 No newline at end of file
+0 −182
Original line number Diff line number Diff line
import math
import unittest
import logging
import re

import numpy as np
from imageio import imread
from scipy.ndimage.interpolation import rotate

from autocnet.examples import get_path
from autocnet.transformation import roi
from .. import ciratefi

import pytest

# Can be parameterized for more exhaustive tests
upsampling = 10
alpha = math.pi/2
cifi_thresh = 90
rafi_thresh = 90
tefi_thresh = 100
use_percentile = True
radii = list(range(1, 3))

@pytest.fixture
def img():
    return imread(get_path('AS15-M-0298_SML.png'), as_gray=True)

@pytest.fixture
def img_coord():
    return 482.09783936, 652.40679932

@pytest.fixture
def template(img, img_coord):
    coord_x, coord_y = img_coord
    template= roi.Roi(img, coord_x, coord_y, 5, 5).clip()
    template = rotate(template, 90)
    return template

@pytest.fixture
def search():
    coord_x, coord_y = (482.09783936, 652.40679932)
    img = imread(get_path('AS15-M-0298_SML.png'), as_gray=True)
    search = roi.Roi(img, coord_x, coord_y, 21, 21).clip()
    return search

@pytest.fixture
def offset_template(img, img_coord):
    coord_x, coord_y = img_coord
    coord_x += 1
    coord_y += 1
    offset_template = roi.Roi(img, coord_x, coord_y, 5, 5).clip()
    return offset_template

def test_cifi_radii_too_large(template, search, caplog):
    # check all logs
    ciratefi.cifi(template, search, 1.0, radii=[100], use_percentile=False)
    num_pattern = '\d+'
    captured_log = caplog.records[0].getMessage()
    match = re.findall(num_pattern, captured_log)
    assert (f'Max Radii is larger than original template, '\
        f'this may produce sub-par results.Max radii: {match[0]} max template dimension: {match[1]}'\
             == caplog.records[0].getMessage())


def test_cifi_bounds_error(template, search):
    with pytest.raises(ValueError), pytest.warns(UserWarning):
        ciratefi.cifi(template, search, -1.1, use_percentile=False)

def test_cifi_radii_none_error(template, search):
    with pytest.raises(ValueError):
        ciratefi.cifi(template, search, 90, radii=None)

def test_cifi_scales_none_error(template, search):
    with pytest.raises(ValueError):
        ciratefi.cifi(template, search, 90, scales=None)

def test_cifi_template_too_large_error(template, search):
    with pytest.raises(ValueError):
        ciratefi.cifi(search,template, 90, scales=None)

@pytest.mark.parametrize('cifi_thresh, radii', [(90,list(range(1, 3)))])
@pytest.mark.filterwarnings('ignore::UserWarning')  # skimage deprecation warnings to move to new defaults
def test_cifi(template, search, cifi_thresh, radii):
    pixels, scales = ciratefi.cifi(template, search, thresh=cifi_thresh,
                                   radii=radii, use_percentile=True)

    assert search.shape == scales.shape
    assert (np.floor(search.shape[0]/2), np.floor(search.shape[1]/2)) in pixels
    assert pixels.size in range(0,search.size)

@pytest.mark.filterwarnings('ignore::UserWarning')  # skimage deprecation warnings to move to new defaults
def test_rafi_warning(template, search,caplog):
    rafi_pixels = [(10, 10)]
    rafi_scales = np.ones(search.shape, dtype=float)
    ciratefi.rafi(template, search, rafi_pixels,
            rafi_scales, thresh=1, radii=[100],
            use_percentile=False)

def test_rafi_bounds_error(template, search):
    rafi_pixels = [(10, 10)]
    rafi_scales = np.ones(search.shape, dtype=float)
    with pytest.raises(ValueError) as f, pytest.warns(UserWarning) as g:
        ciratefi.rafi(template, search, rafi_pixels, rafi_scales, -1.1, use_percentile=False)

def test_rafi_radii_list_none_error(template, search):
    rafi_pixels = [(10, 10)]
    with pytest.raises(ValueError):
        ciratefi.rafi(search, template, rafi_pixels, -1.1, radii=None)

def test_rafi_pixel_list_error(template, search):
    rafi_pixels = []
    rafi_scales = np.ones(search.shape, dtype=float)
    with pytest.raises(ValueError):
        ciratefi.rafi(template, search, rafi_pixels, rafi_scales)

def test_rafi_scales_list_error(template, search):
    rafi_pixels = [(10, 10)]
    with pytest.raises(ValueError):
        ciratefi.rafi(template, search, rafi_pixels, None)

def test_rafi_template_bigger_error(template, search):
    rafi_pixels = [(10, 10)]
    rafi_scales = np.ones(search.shape, dtype=float)
    with pytest.raises(ValueError):
        ciratefi.rafi(search, template, rafi_pixels,rafi_scales)

def test_rafi_shape_mismatch(template, search):
        rafi_pixels = [(10, 10)]
        rafi_scales = np.ones(search.shape, dtype=float)[:10]
        with pytest.raises(ValueError):
            ciratefi.rafi(template, search, rafi_pixels, rafi_scales)

@pytest.mark.parametrize("rafi_thresh, radii, alpha", [(90, list(range(1, 3)),math.pi/2)])
@pytest.mark.filterwarnings('ignore::UserWarning')  # skimage deprecation warnings to move to new defaults
def test_rafi(template, search, rafi_thresh, radii, alpha):
    rafi_pixels = [(10, 10)]
    rafi_scales = np.ones(search.shape, dtype=float)
    pixels, scales = ciratefi.rafi(template, search, rafi_pixels, rafi_scales,
                                   thresh=rafi_thresh, radii=radii, use_percentile=True,
                                   alpha=alpha)

    assert (np.floor(search.shape[0]/4), np.floor(search.shape[1]/4)) in pixels
    assert pixels.size in range(0, search.size)

# Alternate approach to the more verbose tests above - this tests all combinations
@pytest.mark.parametrize("tefi_pixels", [[(10,10)], None])
@pytest.mark.parametrize("tefi_scales", [np.ones((42,42), dtype=float),
                                         None,
                                         np.ones((42,42), dtype=float)[:10]])
@pytest.mark.parametrize("tefi_angles", [[3.14159265], None])
@pytest.mark.parametrize("thresh", [-1.1, 90])
@pytest.mark.parametrize("reverse", [False, True])
def test_tefi_errors(template, search, tefi_pixels, tefi_scales, tefi_angles, thresh, reverse):
    with pytest.raises(ValueError):
        if reverse:
            template, search = search, template
        ciratefi.tefi(template, search, tefi_pixels, tefi_scales,
                  tefi_angles, thresh=-1.1, use_percentile=False, alpha=math.pi/2)

@pytest.mark.filterwarnings('ignore::UserWarning')  # skimage deprecation warnings to move to new defaults
def test_tefi(template, search):
    tefi_pixels = [(10, 10)]
    tefi_scales = np.ones(search.shape, dtype=float)
    tefi_angles = [3.14159265]

    pixel = ciratefi.tefi(template, search, tefi_pixels, tefi_scales, tefi_angles,
                                   thresh=tefi_thresh, use_percentile=True, alpha=math.pi/2,
                                   upsampling=10)

    assert np.equal((11.5, 11.5), (pixel[1], pixel[0])).all()

@pytest.mark.filterwarnings('ignore::UserWarning')  # skimage deprecation warnings to move to new defaults
@pytest.mark.parametrize("cifi_thresh, rafi_thresh, tefi_thresh, alpha, radii",[
                         (90,90,100,math.pi/2,list(range(1, 3)))])
def test_ciratefi(template, search, cifi_thresh, rafi_thresh, tefi_thresh, alpha, radii):
    results = ciratefi.ciratefi(template, search, upsampling=10, cifi_thresh=cifi_thresh,
                                rafi_thresh=rafi_thresh, tefi_thresh=tefi_thresh,
                                use_percentile=True, alpha=alpha, radii=radii)

    assert len(results) == 3
    assert (np.array(results[1], results[0]) < 1).all()
+0 −12
Original line number Diff line number Diff line
@@ -20,15 +20,3 @@ def test_bad_mi():
    corrilation = mutual_information.mutual_information(test_image1, test_image2)
    assert corrilation == pytest.approx(0)
def test_mutual_information():
    d_template = np.array([[i for i in range(50, 100)] for j in range(50)])
    s_image = np.ones((100, 100))

    s_image[25:75, 25:75] = d_template

    x_offset, y_offset, max_corr, corr_map = mutual_information.mutual_information_match(d_template, s_image, bins=20)
    assert x_offset == 0.01711861257171421
    assert y_offset == 0.0
    assert max_corr == 2.9755967600033015
    assert corr_map.shape == (51, 51)
    assert np.min(corr_map) >= 0.0
Loading