From 11323a15279d6b3936c9158aa17c3f868253553b Mon Sep 17 00:00:00 2001 From: Nathan Faggian Date: Sat, 17 Mar 2012 20:27:55 +1100 Subject: [PATCH 1/8] Update README.rst --- README.rst | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/README.rst b/README.rst index 633a9b0..5a3e7ac 100644 --- a/README.rst +++ b/README.rst @@ -3,14 +3,14 @@ About ===== -python-register is a python module for image registration built ontop of scipy and numpy. +imreg, short for image registration, is a python package for image registration built ontop of scipy and numpy. It is currently maintained by Nathan Faggian, Riaan Van Den Dool and Stefan Van Der Walt. Important links =============== -- Official source code: https://github.com/nfaggian/python-register +- Forked from : https://github.com/nfaggian/python-register Dependencies ============ @@ -20,7 +20,6 @@ setuptools, NumPy >= 1.5, SciPy >= 0.9 and a working C++ compiler. To run the tests you will also need py.test >= 2.0. - Install ======= @@ -38,7 +37,7 @@ To install for all users on Unix/Linux:: Development =========== -Basic rules for commits to the python-register repository: +Basic rules for commits to the imreg repository: + master is our stable "release" branch. @@ -51,7 +50,7 @@ GIT You can check the latest sources with the command:: - git clone git://github.com/nfaggian/python-regsiter.git + git clone git://github.com/pyimreg/imreg.git Contributors ~~~~~~~~~~~~~ From ddf8fc78cffbeb809551003c37ba4c61fcba236d Mon Sep 17 00:00:00 2001 From: Nathan Faggian Date: Sat, 17 Mar 2012 20:28:28 +1100 Subject: [PATCH 2/8] Update README.rst --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 5a3e7ac..5b25176 100644 --- a/README.rst +++ b/README.rst @@ -3,7 +3,7 @@ About ===== -imreg, short for image registration, is a python package for image registration built ontop of scipy and numpy. +*imreg*, short for image registration, is a python package for image registration built ontop of scipy and numpy. It is currently maintained by Nathan Faggian, Riaan Van Den Dool and Stefan Van Der Walt. From ee010cee3d68e976b5e66ce1f9164e7263f850b0 Mon Sep 17 00:00:00 2001 From: Nathan Faggian Date: Sat, 17 Mar 2012 20:42:13 +1100 Subject: [PATCH 3/8] Update LICENSE.txt --- LICENSE.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/LICENSE.txt b/LICENSE.txt index 08218b3..856d9e6 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,4 +1,5 @@ -Copyright 2011, Nathan Faggian (nathan.faggian@gmail.com) +Copyright 2011, Nathan Faggian, Austratlian Bureau of Meteorology. +(nfaggian@bom.gov.au) Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. From 40d978397ff2728fcc3588d3e8c5f9f16f583ec9 Mon Sep 17 00:00:00 2001 From: Nathan Faggian Date: Sat, 17 Mar 2012 20:44:09 +1100 Subject: [PATCH 4/8] Update CONTRIBUTORS.txt --- CONTRIBUTORS.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index b17c112..f9654d8 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -1,8 +1,8 @@ - Nathan Faggian - BDFL: Project coordination, methods development. + Australian Bureau of Meteorology - Stefan Van Der Walt - Project coordination, methods development. + UC Berkley - Riaan Van Den Dool - Feature detection methods. \ No newline at end of file + \ No newline at end of file From b57e9e7029ccf8547f905cf2fa2dc75ee7485f24 Mon Sep 17 00:00:00 2001 From: Nathan Faggian Date: Mon, 19 Mar 2012 15:37:52 +1100 Subject: [PATCH 5/8] Update LICENSE.txt --- LICENSE.txt | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/LICENSE.txt b/LICENSE.txt index 856d9e6..384e40c 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,5 +1,8 @@ -Copyright 2011, Nathan Faggian, Austratlian Bureau of Meteorology. -(nfaggian@bom.gov.au) +Copyright 2011, +Nathan Faggian, Australian Bureau of Meteorology. (nfaggian@bom.gov.au) +Stefan Van Der Walt. +Riaan Van Den Dool. + Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. From 4bfe20dc8dfd80fff2f03c19aed3e831e739a0c5 Mon Sep 17 00:00:00 2001 From: Riaan van den Dool Date: Mon, 19 Mar 2012 08:52:28 +0200 Subject: [PATCH 6/8] Added my email address. Do not see a reason to include affiliation info. --- LICENSE.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LICENSE.txt b/LICENSE.txt index 384e40c..f4f65c3 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,7 +1,7 @@ Copyright 2011, Nathan Faggian, Australian Bureau of Meteorology. (nfaggian@bom.gov.au) Stefan Van Der Walt. -Riaan Van Den Dool. +Riaan Van Den Dool. (riaanvddool@gmail.com) Licensed under the Apache License, Version 2.0 (the "License"); From 0a9ea732cbb075865e78145b109bcc3c143386c9 Mon Sep 17 00:00:00 2001 From: Riaan van den Dool Date: Mon, 19 Mar 2012 15:45:34 +0200 Subject: [PATCH 7/8] Corrected some mistakes for documentation purpose. --- documentation/deformation_models.rst | 2 +- register/models/model.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/documentation/deformation_models.rst b/documentation/deformation_models.rst index e998aea..a778e85 100644 --- a/documentation/deformation_models.rst +++ b/documentation/deformation_models.rst @@ -2,5 +2,5 @@ Deformation models ============================= -.. automodule:: register +.. automodule:: register.models.model :members: diff --git a/register/models/model.py b/register/models/model.py index ee8e683..059f6d8 100644 --- a/register/models/model.py +++ b/register/models/model.py @@ -10,15 +10,15 @@ class Model(object): Attributes ---------- - METRIC : string - The type of similarity metric being used. + MODEL : string + The type of deformation model being used. DESCRIPTION : string A meaningful description of the model used, with references where appropriate. """ - MODEL=None - DESCRIPTION=None + MODEL='' + DESCRIPTION='' def __init__(self, coordinates): From 6ab79fa2edfbf0629f17247ef3046f59cb9fe31d Mon Sep 17 00:00:00 2001 From: Nathan Faggian Date: Thu, 22 Mar 2012 17:09:58 +1100 Subject: [PATCH 8/8] First step in refactor. --- documentation/conf.py | 10 +- examples/README.rst | 2 +- examples/haar_features.py | 2 +- examples/linreg.py | 10 +- examples/linreg_pyramid.py | 10 +- examples/nonlinreg.py | 10 +- examples/nonlinreg_image.py | 10 +- examples/nonlinreg_image_features.py | 8 +- examples/nonlinreg_pyramid.py | 10 +- examples/projective.py | 10 +- register/__init__.py | 0 register/features/__init__.py | 0 register/features/detector.py | 96 -- register/features/haar2d.py | 128 --- register/metrics/__init__.py | 0 register/metrics/metric.py | 135 --- register/models/.gitignore | 0 register/models/__init__.py | 0 register/models/model.py | 921 ------------------ register/register.py | 593 ----------- register/samplers/__init__.py | 0 .../lib.macosx-10.6-x86_64-2.7/libsampler.so | Bin 9528 -> 0 bytes .../temp.macosx-10.6-x86_64-2.7/libsampler.o | Bin 48620 -> 0 bytes register/samplers/libsampler.cpp | 285 ------ register/samplers/ndarray.h | 206 ---- register/samplers/numpyctypes.py | 109 --- register/samplers/sampler.py | 266 ----- register/samplers/setup.py | 24 - register/setup.py | 21 - register/visualize/__init__.py | 0 register/visualize/dialog.py | 158 --- register/visualize/dialog.ui | 80 -- register/visualize/plot.py | 225 ----- setup.py | 7 +- tests/test_detectors.py | 2 +- tests/test_haar2d.py | 2 +- tests/test_models.py | 6 +- tests/test_register.py | 8 +- tests/test_register_data.py | 2 +- tests/test_register_features.py | 6 +- tests/test_samplers.py | 6 +- 41 files changed, 60 insertions(+), 3308 deletions(-) delete mode 100644 register/__init__.py delete mode 100644 register/features/__init__.py delete mode 100644 register/features/detector.py delete mode 100644 register/features/haar2d.py delete mode 100644 register/metrics/__init__.py delete mode 100644 register/metrics/metric.py delete mode 100644 register/models/.gitignore delete mode 100644 register/models/__init__.py delete mode 100644 register/models/model.py delete mode 100644 register/register.py delete mode 100644 register/samplers/__init__.py delete mode 100755 register/samplers/build/lib.macosx-10.6-x86_64-2.7/libsampler.so delete mode 100644 register/samplers/build/temp.macosx-10.6-x86_64-2.7/libsampler.o delete mode 100644 register/samplers/libsampler.cpp delete mode 100644 register/samplers/ndarray.h delete mode 100644 register/samplers/numpyctypes.py delete mode 100644 register/samplers/sampler.py delete mode 100644 register/samplers/setup.py delete mode 100644 register/setup.py delete mode 100644 register/visualize/__init__.py delete mode 100644 register/visualize/dialog.py delete mode 100644 register/visualize/dialog.ui delete mode 100644 register/visualize/plot.py diff --git a/documentation/conf.py b/documentation/conf.py index cf31a96..b682c2c 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -40,8 +40,8 @@ master_doc = 'index' # General information about the project. -project = u'python-register' -copyright = u'2011, Nathan Faggian, Stefan Van Der Walt, Riaan Van Den Dool' +project = u'imreg' +copyright = u'2012, Nathan Faggian, Stefan Van Der Walt, Riaan Van Den Dool' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -164,7 +164,7 @@ #html_file_suffix = None # Output file base name for HTML help builder. -htmlhelp_basename = 'python-registerdoc' +htmlhelp_basename = 'imregdoc' # -- Options for LaTeX output -------------------------------------------------- @@ -178,7 +178,7 @@ # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). latex_documents = [ - ('index', 'python-register.tex', u'python-register Documentation', + ('index', 'imreg.tex', u'imreg Documentation', u'Nathan Faggian, Stefan Van Der Walt, Riaan Van Den Dool', 'manual'), ] @@ -211,6 +211,6 @@ # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). man_pages = [ - ('index', 'python-register', u'python-register Documentation', + ('index', 'imreg', u'python-register Documentation', [u'Nathan Faggian, Stefan Van Der Walt, Riaan Van Den Dool'], 1) ] diff --git a/examples/README.rst b/examples/README.rst index da591f2..1c5bfa6 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -3,7 +3,7 @@ About ======== -Demonstrates uses of the sci-kit: +Demonstrates uses of the package: Examples of 2D linear and nonlinear image registration. diff --git a/examples/haar_features.py b/examples/haar_features.py index f4ef46f..f42d2dc 100644 --- a/examples/haar_features.py +++ b/examples/haar_features.py @@ -6,7 +6,7 @@ import scipy.ndimage as nd from matplotlib.pyplot import imread, plot, imshow, show -from register.features.detector import detect, HaarDetector +from imreg.features.detector import detect, HaarDetector # Load the image. image = imread('data/cameraman.png') diff --git a/examples/linreg.py b/examples/linreg.py index cd68a60..1a398cb 100644 --- a/examples/linreg.py +++ b/examples/linreg.py @@ -7,12 +7,12 @@ import scipy.ndimage as nd import scipy.misc as misc -from register.models import model -from register.metrics import metric -from register.samplers import sampler -from register import register +from imreg.models import model +from imreg.metrics import metric +from imreg.samplers import sampler +from imreg import register -from register.visualize import plot +from imreg.visualize import plot # Form some test data (lena, lena rotated 20 degrees) image = misc.lena() diff --git a/examples/linreg_pyramid.py b/examples/linreg_pyramid.py index 1bcb357..6ad39b9 100644 --- a/examples/linreg_pyramid.py +++ b/examples/linreg_pyramid.py @@ -6,13 +6,13 @@ import scipy.ndimage as nd import scipy.misc as misc -from register.models import model -from register.metrics import metric -from register.samplers import sampler +from imreg.models import model +from imreg.metrics import metric +from imreg.samplers import sampler -from register.visualize import plot +from imreg.visualize import plot -from register import register +from imreg import register # Form some test data (lena, lena rotated 20 degrees) image = register.RegisterData(misc.lena()) diff --git a/examples/nonlinreg.py b/examples/nonlinreg.py index 44b8903..6951a71 100644 --- a/examples/nonlinreg.py +++ b/examples/nonlinreg.py @@ -7,12 +7,12 @@ import scipy.ndimage as nd import scipy.misc as misc -from register.models import model -from register.metrics import metric -from register.samplers import sampler +from imreg.models import model +from imreg.metrics import metric +from imreg.samplers import sampler -from register.visualize import plot -from register import register +from imreg.visualize import plot +from imreg import register def warp(image): """ diff --git a/examples/nonlinreg_image.py b/examples/nonlinreg_image.py index 8d19e99..12f20fa 100644 --- a/examples/nonlinreg_image.py +++ b/examples/nonlinreg_image.py @@ -11,12 +11,12 @@ from matplotlib.pyplot import imread -from register.models import model -from register.metrics import metric -from register.samplers import sampler +from imreg.models import model +from imreg.metrics import metric +from imreg.samplers import sampler -from register.visualize import plot -from register import register +from imreg.visualize import plot +from imreg import register # Form some test data (lena, lena rotated 20 degrees) image = imread('data/frown.png')[:, :, 0] diff --git a/examples/nonlinreg_image_features.py b/examples/nonlinreg_image_features.py index 61577d7..e66d9ff 100644 --- a/examples/nonlinreg_image_features.py +++ b/examples/nonlinreg_image_features.py @@ -13,11 +13,11 @@ import matplotlib.pyplot as plt -from register.models import model -from register.samplers import sampler -from register.visualize import plot +from imreg.models import model +from imreg.samplers import sampler +from imreg.visualize import plot -from register import register +from imreg import register # Load the image and feature data. image = register.RegisterData( diff --git a/examples/nonlinreg_pyramid.py b/examples/nonlinreg_pyramid.py index 69be0b4..1c60bd7 100644 --- a/examples/nonlinreg_pyramid.py +++ b/examples/nonlinreg_pyramid.py @@ -7,11 +7,11 @@ import scipy.ndimage as nd import scipy.misc as misc -from register.models import model -from register.metrics import metric -from register.samplers import sampler -from register.visualize import plot -from register import register +from imreg.models import model +from imreg.metrics import metric +from imreg.samplers import sampler +from imreg.visualize import plot +from imreg import register def warp(image): """ diff --git a/examples/projective.py b/examples/projective.py index 9135f91..0b27b98 100644 --- a/examples/projective.py +++ b/examples/projective.py @@ -7,12 +7,12 @@ import scipy.ndimage as nd import scipy.misc as misc -from register.models import model -from register.metrics import metric -from register.samplers import sampler +from imreg.models import model +from imreg.metrics import metric +from imreg.samplers import sampler -from register.visualize import plot -from register import register +from imreg.visualize import plot +from imreg import register def warp(image): """ diff --git a/register/__init__.py b/register/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/features/__init__.py b/register/features/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/features/detector.py b/register/features/detector.py deleted file mode 100644 index 6d2fabc..0000000 --- a/register/features/detector.py +++ /dev/null @@ -1,96 +0,0 @@ -""" A collection of feature detectors""" - -import numpy as np -import scipy.ndimage as nd -import math - -from register.features.haar2d import haar2d - -__debug=False - - -def _debug(something): - global __debug - if __debug: - print something - - -# Constants -HaarDetector = 0 - -def detect(image, detectorType=HaarDetector, options=None, debug=False): - global __debug - global __plt - __debug = debug - if detectorType == HaarDetector: - return _detectHaarFeatures(image, options) - else: # default detector - return _detectHaarFeatures(image, options) - - -def _haarDefaultOptions(image): - options = {} - options['levels'] = 5 # number of wavelet levels - options['threshold'] = 0.2 # threshold between 0.0 and 1.0 to filter out weak features (0.0 includes all features) - options['locality'] = 5 # minimum (approx) distance between two features - return options - -def _detectHaarFeatures(image, options={}): - if options is None: - options = _haarDefaultOptions(image) - levels = options.get('levels') - maxpoints = options.get('maxpoints') - threshold = options.get('threshold') - locality = options.get('locality') - - haarData = haar2d(image, levels) - - avgRows = haarData.shape[0] / 2 ** levels - avgCols = haarData.shape[1] / 2 ** levels - - SalientPoints = {} - - siloH = np.zeros([haarData.shape[0]/2, haarData.shape[1]/2, levels]) - siloD = np.zeros([haarData.shape[0]/2, haarData.shape[1]/2, levels]) - siloV = np.zeros([haarData.shape[0]/2, haarData.shape[1]/2, levels]) - - # Build the saliency silos - for i in range(levels): - level = i + 1 - halfRows = haarData.shape[0] / 2 ** level - halfCols = haarData.shape[1] / 2 ** level - siloH[:,:,i] = nd.zoom(haarData[:halfRows, halfCols:halfCols*2], 2**(level-1)) - siloD[:,:,i] = nd.zoom(haarData[halfRows:halfRows*2, halfCols:halfCols*2], 2**(level-1)) - siloV[:,:,i] = nd.zoom(haarData[halfRows:halfRows*2, :halfCols], 2**(level-1)) - - # Calculate saliency heat-map - saliencyMap = np.max(np.array([ - np.sum(np.abs(siloH), axis=2), - np.sum(np.abs(siloD), axis=2), - np.sum(np.abs(siloV), axis=2) - ]), axis=0) - - # Determine global maximum and saliency threshold - maximum = np.max(saliencyMap) - sthreshold = threshold * maximum - - # Extract features by finding local maxima - rows = haarData.shape[0] / 2 - cols = haarData.shape[1] / 2 - features = {} - id = 0 - for row in range(locality,rows-locality): - for col in range(locality,cols-locality): - saliency = saliencyMap[row,col] - if saliency > sthreshold: - if saliency >= np.max(saliencyMap[row-locality:row+locality, col-locality:col+locality]): - features[id] = (row*2,col*2) - id += 1 - - result = {} - result['points'] = features - return result - - - - diff --git a/register/features/haar2d.py b/register/features/haar2d.py deleted file mode 100644 index 5b88c0c..0000000 --- a/register/features/haar2d.py +++ /dev/null @@ -1,128 +0,0 @@ -import numpy as np -import scipy.ndimage as nd -import math - -__debug=False - -def _debug(something): - global __debug - if __debug: - print something - -def haar2d(image, levels, debug=False): - """ - 2D Haar wavelet decomposition for levels=levels. - - Parameters - ---------- - image: nd-array - Input image - levels: int - Number of wavelet levels to compute - debug: - Setting debug=True will produce some debug output messages - - Returns - ------- - haarImage: nd-array - An image containing the Haar decomposition of the input image. - Might be larger than the input image. - - See also - -------- - register.features.ihaar2d - """ - global __debug - __debug = debug - assert len(image.shape) == 2, 'Must be 2D image!' - origRows, origCols = image.shape - extraRows = 0; - extraCols = 0; - while (((origRows + extraRows) >> levels) << levels != (origRows + extraRows)): - extraRows += 1 - while (((origCols + extraCols) >> levels) << levels != (origCols + extraCols)): - extraCols += 1 - _debug("Padding: %d x %d -> %d x %d" % (origRows, origCols, origRows + extraRows, origCols + extraCols)) - - # Pad image to compatible shape using repitition - rightFill = np.repeat(image[:, -1:], extraCols, axis=1) - _image = np.zeros([origRows, origCols + extraCols]) - _image[:, :origCols] = image - _image[:, origCols:] = rightFill - bottomFill = np.repeat(_image[-1:, :], extraRows, axis=0) - image = np.zeros([origRows + extraRows, origCols + extraCols]) - image[:origRows, :] = _image - image[origRows:, :] = bottomFill - _debug("Padded image is: %d x %d" % (image.shape[0], image.shape[1])) - - haarImage = image - for level in range(1,levels+1): - halfRows = image.shape[0] / 2 ** level - halfCols = image.shape[1] / 2 ** level - _image = image[:halfRows*2, :halfCols*2] - # rows - lowpass = (_image[:, :-1:2] + _image[:, 1::2]) / 2 - higpass = (_image[:, :-1:2] - _image[:, 1::2]) / 2 - _image[:, :_image.shape[1]/2] = lowpass - _image[:, _image.shape[1]/2:] = higpass - # cols - lowpass = (_image[:-1:2, :] + _image[1::2, :]) / 2 - higpass = (_image[:-1:2, :] - _image[1::2, :]) / 2 - _image[:_image.shape[0]/2, :] = lowpass - _image[_image.shape[0]/2:, :] = higpass - haarImage[:halfRows*2, :halfCols*2] = _image - - _debug(haarImage) - return haarImage - -def ihaar2d(image, levels, debug=False): - """ - 2D Haar wavelet decomposition inverse for levels=levels. - - Parameters - ---------- - image: nd-array - Input image - levels: int - Number of wavelet levels to de-compute - debug: - Setting debug=True will produce some debug output messages - - Returns - ------- - image: nd-array - An image containing the inverse Haar decomposition of the input image. - - See also - -------- - register.features.haar2d - """ - global __debug - __debug = debug - assert len(image.shape) == 2, 'Must be 2D image!' - origRows, origCols = image.shape - extraRows = 0; - extraCols = 0; - while (((origRows + extraRows) >> levels) << levels != (origRows + extraRows)): - extraRows += 1 - while (((origCols + extraCols) >> levels) << levels != (origCols + extraCols)): - extraCols += 1 - assert (extraRows, extraCols) == (0,0), 'Must be compatible shape!' - - for level in range(levels, 0, -1): - _debug("level=%d" % level) - halfRows = image.shape[0] / 2 ** level - halfCols = image.shape[1] / 2 ** level - # cols - lowpass = image[:halfRows*2, :halfCols].copy() - higpass = image[:halfRows*2, halfCols:halfCols*2].copy() - image[:halfRows*2, :halfCols*2-1:2] = lowpass + higpass - image[:halfRows*2, 1:halfCols*2:2] = lowpass - higpass - _debug(image) - # rows - lowpass = image[:halfRows, :halfCols*2].copy() - higpass = image[halfRows:halfRows*2, :halfCols*2].copy() - image[:halfRows*2-1:2, :halfCols*2] = lowpass + higpass - image[1:halfRows*2:2, :halfCols*2] = lowpass - higpass - - return image diff --git a/register/metrics/__init__.py b/register/metrics/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/metrics/metric.py b/register/metrics/metric.py deleted file mode 100644 index 5bce219..0000000 --- a/register/metrics/metric.py +++ /dev/null @@ -1,135 +0,0 @@ -""" A collection of image similarity metrics. """ - -import numpy as np - -class Metric(object): - """ - Abstract similarity metric. - - Attributes - ---------- - METRIC : string - The type of similarity metric being used. - DESCRIPTION : string - A meaningful description of the metric used, with references where - appropriate. - """ - - METRIC=None - DESCRIPTION=None - - def __init__(self): - pass - - def error(self, warpedImage, template): - """ - Evaluates the metric. - - Parameters - ---------- - warpedImage: nd-array - Input image after warping. - template: nd-array - Template image. - - Returns - ------- - error: nd-array - Metric evaluated over all image coordinates. - """ - - raise NotImplementedError('') - - def jacobian(self, model, warpedImage, p=None): - """ - Computes the jacobian dP/dE. - - Parameters - ---------- - model: deformation model - A particular deformation model. - warpedImage: nd-array - Input image after warping. - p : optional list - Current warp parameters - - Returns - ------- - jacobian: nd-array - A derivative of model parameters with respect to the metric. - """ - raise NotImplementedError('') - - def __str__(self): - return 'Metric: {0} \n {1}'.format( - self.METRIC, - self.DESCRIPTION - ) - - -class Residual(Metric): - """ Standard least squares metric """ - - METRIC='residual' - - DESCRIPTION=""" - The residual which is computed as the difference between the - deformed image an the template: - - (I(W(x;p)) - T) - - """ - - def __init__(self): - Metric.__init__(self) - - def jacobian(self, model, warpedImage, p=None): - """ - Computes the jacobian dP/dE. - - Parameters - ---------- - model: deformation model - A particular deformation model. - warpedImage: nd-array - Input image after warping. - p : optional list - Current warp parameters - - Returns - ------- - jacobian: nd-array - A jacobain matrix. (m x n) - | where: m = number of image pixels, - | p = number of parameters. - """ - - grad = np.gradient(warpedImage) - - dIx = grad[1].flatten() - dIy = grad[0].flatten() - - dPx, dPy = model.jacobian(p) - - J = np.zeros_like(dPx) - for index in range(0, dPx.shape[1]): - J[:,index] = dPx[:,index]*dIx + dPy[:,index]*dIy - return J - - def error(self, warpedImage, template): - """ - Evaluates the residual metric. - - Parameters - ---------- - warpedImage: nd-array - Input image after warping. - template: nd-array - Template image. - - Returns - ------- - error: nd-array - Metric evaluated over all image coordinates. - """ - return warpedImage.flatten() - template.flatten() diff --git a/register/models/.gitignore b/register/models/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/register/models/__init__.py b/register/models/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/models/model.py b/register/models/model.py deleted file mode 100644 index 059f6d8..0000000 --- a/register/models/model.py +++ /dev/null @@ -1,921 +0,0 @@ -""" A collection of deformation models. """ - -import numpy as np -import scipy.signal as signal - - -class Model(object): - """ - Abstract geometry model. - - Attributes - ---------- - MODEL : string - The type of deformation model being used. - DESCRIPTION : string - A meaningful description of the model used, with references where - appropriate. - """ - - MODEL='' - DESCRIPTION='' - - - def __init__(self, coordinates): - self.coordinates = coordinates - - - def fit(self, p0, p1): - """ - Estimates the best fit parameters that define a warp field, which - deforms feature points p0 to p1. - - Parameters - ---------- - p0: nd-array - Image features (points). - p1: nd-array - Template features (points). - - Returns - ------- - parameters: nd-array - Model parameters. - error: float - Sum of RMS error between p1 and alinged p0. - """ - raise NotImplementedError('') - - - @staticmethod - def scale(p, factor): - """ - Scales an transformtaion by a factor. - - Parameters - ---------- - p: nd-array - Model parameters. - factor: float - A scaling factor. - - Returns - ------- - parameters: nd-array - Model parameters. - """ - raise NotImplementedError('') - - - def estimate(self, warp): - """ - Estimates the best fit parameters that define a warp field. - - Parameters - ---------- - warp: nd-array - Deformation field. - - Returns - ------- - parameters: nd-array - Model parameters. - """ - raise NotImplementedError('') - - - def warp(self, parameters): - """ - Computes the warp field given model parameters. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - warp: nd-array - Deformation field. - """ - - displacement = self.transform(parameters) - - # Approximation of the inverse (samplers work on inverse warps). - return self.coordinates.tensor + displacement - - - def transform(self, parameters): - """ - A geometric transformation of coordinates. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - coords: nd-array - Deformation coordinates. - """ - raise NotImplementedError('') - - - def jacobian(self, p=None): - """ - Evaluates the derivative of deformation model with respect to the - coordinates. - """ - raise NotImplementedError('') - - - def __str__(self): - return 'Model: {0} \n {1}'.format( - self.MODEL, - self.DESCRIPTION - ) - - -class Shift(Model): - - MODEL='Shift (S)' - - DESCRIPTION=""" - Applies the shift coordinate transformation. Follows the derivations - shown in: - - S. Baker and I. Matthews. 2004. Lucas-Kanade 20 Years On: A - Unifying Framework. Int. J. Comput. Vision 56, 3 (February 2004). - """ - - - def __init__(self, coordinates): - Model.__init__(self, coordinates) - - - @property - def identity(self): - return np.zeros(2) - - - @staticmethod - def scale(p, factor): - """ - Scales an shift transformation by a factor. - - Parameters - ---------- - p: nd-array - Model parameters. - factor: float - A scaling factor. - - Returns - ------- - parameters: nd-array - Model parameters. - """ - - pHat = p.copy() - pHat *= factor - return pHat - - - def fit(self, p0, p1, lmatrix=False): - """ - Estimates the best fit parameters that define a warp field, which - deforms feature points p0 to p1. - - Parameters - ---------- - p0: nd-array - Image features (points). - p1: nd-array - Template features (points). - - Returns - ------- - parameters: nd-array - Model parameters. - error: float - Sum of RMS error between p1 and alinged p0. - """ - - parameters = p1.mean(axis=0) - p0.mean(axis=0) - - projP0 = p0 + parameters - - error = np.sqrt( - (projP0[:,0] - p1[:,0])**2 + (projP0[:,1] - p1[:,1])**2 - ).sum() - - return -parameters, error - - - def transform(self, parameters): - """ - A "shift" transformation of coordinates. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - coords: nd-array - Deformation coordinates. - """ - - T = np.eye(3,3) - T[0,2] = -parameters[0] - T[1,2] = -parameters[1] - - displacement = np.dot(T, self.coordinates.homogenous) - \ - self.coordinates.homogenous - - shape = self.coordinates.tensor[0].shape - - return np.array( [ displacement[1].reshape(shape), - displacement[0].reshape(shape) - ] - ) - - - def jacobian(self, p=None): - """ - Evaluates the derivative of deformation model with respect to the - coordinates. - """ - - dx = np.zeros((self.coordinates.tensor[0].size, 2)) - dy = np.zeros((self.coordinates.tensor[0].size, 2)) - - dx[:,0] = 1 - dy[:,1] = 1 - - return (dx, dy) - - -class Affine(Model): - - MODEL='Affine (A)' - - DESCRIPTION=""" - Applies the affine coordinate transformation. Follows the derivations - shown in: - - S. Baker and I. Matthews. 2004. Lucas-Kanade 20 Years On: A - Unifying Framework. Int. J. Comput. Vision 56, 3 (February 2004). - """ - - - def __init__(self, coordinates): - Model.__init__(self, coordinates) - - - @property - def identity(self): - return np.zeros(6) - - - @staticmethod - def scale(p, factor): - """ - Scales an affine transformation by a factor. - - Parameters - ---------- - p: nd-array - Model parameters. - factor: float - A scaling factor. - - Returns - ------- - parameters: nd-array - Model parameters. - """ - - pHat = p.copy() - pHat[4:] *= factor - return pHat - - - def fit(self, p0, p1, lmatrix=False): - """ - Estimates the best fit parameters that define a warp field, which - deforms feature points p0 to p1. - - Parameters - ---------- - p0: nd-array - Image features (points). - p1: nd-array - Template features (points). - - Returns - ------- - parameters: nd-array - Model parameters. - error: float - Sum of RMS error between p1 and alinged p0. - """ - - # Solve: H*X = Y - # --------------------- - # H = Y*inv(X) - - X = np.ones((3, len(p0))) - X[0:2,:] = p0.T - - Y = np.ones((3, len(p0))) - Y[0:2,:] = p1.T - - H = np.dot(Y, np.linalg.pinv(X)) - - parameters = [ - H[0,0] - 1.0, - H[1,0], - H[0,1], - H[1,1] - 1.0, - H[0,2], - H[1,2] - ] - - projP0 = np.dot(H, X)[0:2,:].T - - error = np.sqrt( - (projP0[:,0] - p1[:,0])**2 + (projP0[:,1] - p1[:,1])**2 - ).sum() - - return parameters, error - - - def transform(self, p): - """ - An "affine" transformation of coordinates. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - coords: nd-array - Deformation coordinates. - """ - - T = np.array([ - [p[0]+1.0, p[2], p[4]], - [p[1], p[3]+1.0, p[5]], - [0, 0, 1] - ]) - - displacement = np.dot(np.linalg.inv(T), self.coordinates.homogenous) - \ - self.coordinates.homogenous - - shape = self.coordinates.tensor[0].shape - - return np.array( [ displacement[1].reshape(shape), - displacement[0].reshape(shape) - ] - ) - - - def jacobian(self, p=None): - """" - Evaluates the derivative of deformation model with respect to the - coordinates. - """ - - dx = np.zeros((self.coordinates.tensor[0].size, 6)) - dy = np.zeros((self.coordinates.tensor[0].size, 6)) - - dx[:,0] = self.coordinates.tensor[1].flatten() - dx[:,2] = self.coordinates.tensor[0].flatten() - dx[:,4] = 1.0 - - dy[:,1] = self.coordinates.tensor[1].flatten() - dy[:,3] = self.coordinates.tensor[0].flatten() - dy[:,5] = 1.0 - - return (dx, dy) - - -class Projective(Model): - - MODEL='Projective (P)' - - DESCRIPTION=""" - Applies the projective coordinate transformation. Follows the derivations - shown in: - - S. Baker and I. Matthews. 2004. Lucas-Kanade 20 Years On: A - Unifying Framework. Int. J. Comput. Vision 56, 3 (February 2004). - """ - - - def __init__(self, coordinates): - Model.__init__(self, coordinates) - - - @property - def identity(self): - return np.zeros(9) - - - def fit(self, p0, p1, lmatrix=False): - """ - Estimates the best fit parameters that define a warp field, which - deforms feature points p0 to p1. - - Parameters - ---------- - p0: nd-array - Image features (points). - p1: nd-array - Template features (points). - - Returns - ------- - parameters: nd-array - Model parameters. - error: float - Sum of RMS error between p1 and alinged p0. - """ - - # Solve: H*X = Y - # --------------------- - # H = Y*inv(X) - - X = np.ones((3, len(p0))) - X[0:2,:] = p0.T - - Y = np.ones((3, len(p0))) - Y[0:2,:] = p1.T - - H = np.dot(Y, np.linalg.pinv(X)) - - parameters = [ - H[0,0] - 1.0, - H[1,0], - H[0,1], - H[1,1] - 1.0, - H[0,2], - H[1,2], - H[2,0], - H[2,1], - H[2,2] - 1.0 - ] - - projP0 = np.dot(H, X)[0:2,:].T - - error = np.sqrt( - (projP0[:,0] - p1[:,0])**2 + (projP0[:,1] - p1[:,1])**2 - ).sum() - - return parameters, error - - - def transform(self, p): - """ - An "projective" transformation of coordinates. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - coords: nd-array - Deformation coordinates. - """ - - T = np.array([ - [p[0]+1.0, p[2], p[4]], - [p[1], p[3]+1.0, p[5]], - [p[6], p[7], p[8]+1.0] - ]) - - displacement = np.dot(np.linalg.inv(T), self.coordinates.homogenous) - \ - self.coordinates.homogenous - - shape = self.coordinates.tensor[0].shape - - return np.array( [ displacement[1].reshape(shape), - displacement[0].reshape(shape) - ] - ) - - - def jacobian(self, p): - """" - Evaluates the derivative of deformation model with respect to the - coordinates. - """ - - dx = np.zeros((self.coordinates.tensor[0].size, 9)) - dy = np.zeros((self.coordinates.tensor[0].size, 9)) - - x = self.coordinates.tensor[1].flatten() - y = self.coordinates.tensor[0].flatten() - - dx[:,0] = x / (p[6]*x + p[7]*y + p[8] + 1) - dx[:,2] = y / (p[6]*x + p[7]*y + p[8] + 1) - dx[:,4] = 1.0 / (p[6]*x + p[7]*y + p[8] + 1) - dx[:,6] = x * (p[0]*x + p[2]*y + p[4] + x) / (p[6]*x + p[7]*y + p[8] + 1)**2 - dx[:,7] = y * (p[0]*x + p[2]*y + p[4] + x) / (p[6]*x + p[7]*y + p[8] + 1)**2 - dx[:,8] = 1.0 * (p[0]*x + p[2]*y + p[4] + x) / (p[6]*x + p[7]*y + p[8] + 1)**2 - - dy[:,1] = x / (p[6]*x + p[7]*y + p[8] + 1) - dy[:,3] = y / (p[6]*x + p[7]*y + p[8] + 1) - dy[:,5] = 1.0 / (p[6]*x + p[7]*y + p[8] + 1) - dy[:,6] = x * (p[1]*x + p[3]*y + p[5] + y) / (p[6]*x + p[7]*y + p[8] + 1)**2 - dy[:,7] = y * (p[1]*x + p[3]*y + p[5] + y) / (p[6]*x + p[7]*y + p[8] + 1)**2 - dy[:,8] = 1.0 * (p[1]*x + p[3]*y + p[5] + y) / (p[6]*x + p[7]*y + p[8] + 1)**2 - - return (dx, dy) - - - @staticmethod - def scale(p, factor): - """ - Scales an projective transformation by a factor. - - Derivation: If Hx = x^ , - then SHx = Sx^ , - where S = [[s, 0, 0], [0, s, 0], [0, 0, 1]] . - Now SH = S[[h00, h01, h02], [h10, h11, h12], [h20, h21, h22]] - = [[s.h00, s.h01, s.h02], [s.h10, s.h11, s.h12], [h20, h21, h22]] . - - - Parameters - ---------- - p: nd-array - Model parameters. - factor: float - A scaling factor. - - Returns - ------- - parameters: nd-array - Model parameters. - """ - - pHat = p.copy() - pHat[0:6] *= factor - return pHat - - -class ThinPlateSpline(Model): - - MODEL='Thin Plate Spline (TPS)' - - DESCRIPTION=""" - Computes a thin-plate-spline deformation model, as described in: - - Bookstein, F. L. (1989). Principal warps: thin-plate splines and the - decomposition of deformations. IEEE Transactions on Pattern Analysis - and Machine Intelligence, 11(6), 567-585. - - """ - - def __init__(self, coordinates): - - Model.__init__(self, coordinates) - - def U(self, r): - """ - Kernel function, applied to solve the biharmonic equation. - - Parameters - ---------- - r: float - Distance between sample and coordinate point. - - Returns - ------- - U: float - Evaluated kernel. - """ - - return np.multiply(-np.power(r,2), np.log(np.power(r,2) + 1e-20)) - - def fit(self, p0, p1, lmatrix=False): - """ - Estimates the best fit parameters that define a warp field, which - deforms feature points p0 to p1. - - Parameters - ---------- - p0: nd-array - Image features (points). - p1: nd-array - Template features (points). - lmatrix: boolean - Enables the spline design matrix when returning. - - Returns - ------- - parameters: nd-array - Model parameters. - error: float - Sum of RMS error between p1 and alinged p0. - L: nd-array - Spline design matrix, optional (using lmatrix keyword). - """ - - K = np.zeros((p0.shape[0], p0.shape[0])) - - for i in range(0, p0.shape[0]): - for j in range(0, p0.shape[0]): - r = np.sqrt( (p0[i,0] - p0[j,0])**2 + (p0[i,1] - p0[j,1])**2 ) - K[i,j] = self.U(r) - - P = np.hstack((np.ones((p0.shape[0], 1)), p0)) - - L = np.vstack((np.hstack((K,P)), - np.hstack((P.transpose(), np.zeros((3,3)))))) - - Y = np.vstack( (p1, np.zeros((3, 2))) ) - - parameters = np.dot(np.linalg.inv(L), Y) - - # Estimate the thin-plate spline basis. - self.__basis(p0) - - # Estimate the model fit error. - _p0, _p1, _projP0, error = self.__splineError(p0, p1, parameters) - - if lmatrix: - return parameters, error, L - else: - return parameters, error - - def __splineError(self, p0, p1, parameters): - """ - Estimates the point alignment and computes the alignment error. - - Parameters - ---------- - p0: nd-array - Image features (points). - p1: nd-array - Template features (points). - parameters: nd-array - Thin-plate spline parameters. - - Returns - ------- - error: float - Alignment error between p1 and projected p0 (RMS). - """ - - # like __basis, compute a reduced set of basis vectors. - - basis = np.zeros((p0.shape[0], len(p0)+3)) - - # nonlinear, spline component. - for index, p in enumerate( p0 ): - basis[:,index] = self.U( - np.sqrt( - (p[0]-p1[:,0])**2 + - (p[1]-p1[:,1])**2 - ) - ).flatten() - - # linear, affine component - basis[:,-3] = 1 - basis[:,-2] = p1[:,1] - basis[:,-1] = p1[:,0] - - # compute the alignment error. - - projP0 = np.vstack( [ - np.dot(basis, parameters[:,1]), - np.dot(basis, parameters[:,0]) - ] - ).T - - error = np.sqrt( - (projP0[:,0] - p1[:,0])**2 + (projP0[:,1] - p1[:,1])**2 - ).sum() - - return p0, p1, projP0, error - - def __basis(self, p0): - """ - Forms the thin plate spline deformation basis, which is composed of - a linear and non-linear component. - - Parameters - ---------- - p0: nd-array - Image features (points). - """ - - self.basis = np.zeros((self.coordinates.tensor[0].size, len(p0)+3)) - - # nonlinear, spline component. - for index, p in enumerate( p0 ): - self.basis[:,index] = self.U( - np.sqrt( - (p[0]-self.coordinates.tensor[1])**2 + - (p[1]-self.coordinates.tensor[0])**2 - ) - ).flatten() - - # linear, affine component - - self.basis[:,-3] = 1 - self.basis[:,-2] = self.coordinates.tensor[1].flatten() - self.basis[:,-1] = self.coordinates.tensor[0].flatten() - - - def transform(self, parameters): - """ - A "thin-plate-spline" transformation of coordinates. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - coords: nd-array - Deformation coordinates. - """ - - shape = self.coordinates.tensor[0].shape - - return np.array( [ np.dot(self.basis, parameters[:,1]).reshape(shape), - np.dot(self.basis, parameters[:,0]).reshape(shape) - ] - ) - - def warp(self, parameters): - """ - Computes the warp field given model parameters. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - warp: nd-array - Deformation field. - """ - - return self.transform(parameters) - - - def jacobian(self, p=None): - raise NotImplementedError(""" - It does not make sense to use a non-linear optimization to - fit a thin-plate-spline model. Try the "CubicSpline" deformation - model instead. - """ - ) - - @property - def identity(self): - raise NotImplementedError('') - - -class CubicSpline(Model): - - MODEL='CubicSpline (CS)' - - DESCRIPTION=""" - Applies a cubic-spline deformation model, as described in: - - Kybic, J. and Unser, M. (2003). Fast parametric elastic image - registration. IEEE Transactions on Image Processing, 12(11), 1427-1442. - """ - - def __init__(self, coordinates): - - Model.__init__(self, coordinates) - self.__basis() - - @property - def identity(self): - return np.zeros(self.basis.shape[1]*2) - - @property - def numberOfParameters(self): - return self.basis.shape[1] - - def __basis(self, order=4, divisions=5): - """ - Computes the spline tensor product and stores the products, as basis - vectors. - - Parameters - ---------- - order: int - B-spline order, optional. - divisions: int, optional. - Number of spline knots. - """ - - shape = self.coordinates.tensor[0].shape - grid = self.coordinates.tensor - - spacing = shape[1] / divisions - xKnots = shape[1] / spacing - yKnots = shape[0] / spacing - - Qx = np.zeros((grid[0].size, xKnots)) - Qy = np.zeros((grid[0].size, yKnots)) - - for index in range(0, xKnots): - bx = signal.bspline( grid[1] / spacing - index, order) - Qx[:,index] = bx.flatten() - - for index in range(0, yKnots): - by = signal.bspline( grid[0] / spacing - index, order) - Qy[:,index] = by.flatten() - - basis = [] - for j in range(0,xKnots): - for k in range(0, yKnots): - basis.append(Qx[:,j]*Qy[:,k]) - - self.basis = np.array(basis).T - - - def estimate(self, warp): - """ - Estimates the best fit parameters that define a warp field. - - Parameters - ---------- - warp: nd-array - A deformation field. - - Returns - ------- - parameters: nd-array - Model parameters. - """ - - invB = np.linalg.pinv(self.basis) - - return np.hstack( - (np.dot(invB, warp[1].flatten()), - np.dot(invB, warp[0].flatten())) - ) - - - def transform(self, p): - """ - Applies an spline transformation to image coordinates. - - Parameters - ---------- - parameters: nd-array - Model parameters. - - Returns - ------- - coordinates: nd-array - Deformation coordinates. - """ - - px = np.array(p[0:self.numberOfParameters]) - py = np.array(p[self.numberOfParameters::]) - - shape = self.coordinates.tensor[0].shape - - # FIXME: Inverse of a warp field needs to be derived and put in here, - # clearly a multiplication by -1 is not a good approach. - - return -1.0 * np.array( [ np.dot(self.basis, py).reshape(shape), - np.dot(self.basis, px).reshape(shape) - ] - ) - - def jacobian(self, p=None): - """ - Evaluate the derivative of deformation model with respect to the - coordinates. - """ - - dx = np.zeros((self.coordinates.tensor[0].size, - 2*self.numberOfParameters)) - - dy = np.zeros((self.coordinates.tensor[0].size, - 2*self.numberOfParameters)) - - dx[:, 0:self.numberOfParameters] = self.basis - dy[:, self.numberOfParameters::] = self.basis - - return (dx, dy) diff --git a/register/register.py b/register/register.py deleted file mode 100644 index f10f52c..0000000 --- a/register/register.py +++ /dev/null @@ -1,593 +0,0 @@ -""" A top level registration module """ - -import numpy as np -import scipy.ndimage as nd - -def _smooth(image, variance): - """ - Gaussian smoothing using the fast-fourier-transform (FFT) - - Parameters - ---------- - image: nd-array - Input image - variance: float - Variance of the Gaussian kernel. - - Returns - ------- - image: nd-array - An image convolved with the Gaussian kernel. - - See also - -------- - regisger.Register.smooth - """ - - return np.real( - np.fft.ifft2( - nd.fourier_gaussian( - np.fft.fft2(image), - variance - ) - ) - ) - - -class Coordinates(object): - """ - Container for grid coordinates. - - Attributes - ---------- - domain : nd-array - Domain of the coordinate system. - tensor : nd-array - Grid coordinates. - homogenous : nd-array - `Homogenous` coordinate system representation of grid coordinates. - """ - - def __init__(self, domain, spacing=None): - - self.spacing = 1.0 if not spacing else spacing - self.domain = domain - self.tensor = np.mgrid[0.:domain[1], 0.:domain[3]] - - self.homogenous = np.zeros((3,self.tensor[0].size)) - self.homogenous[0] = self.tensor[1].flatten() - self.homogenous[1] = self.tensor[0].flatten() - self.homogenous[2] = 1.0 - - -class RegisterData(object): - """ - Container for registration data. - - Attributes - ---------- - data : nd-array - The image registration image values. - coords : nd-array, optional - The grid coordinates. - features : dictionary, optional - A mapping of unique ids to registration features. - """ - - def __init__(self, data, coords=None, features=None, spacing=1.0): - - self.data = data.astype(np.double) - - if not coords: - self.coords = Coordinates( - [0, data.shape[0], 0, data.shape[1]], - spacing=spacing - ) - else: - self.coords = coords - - # Features are (as a starting point a dictionary) which define - # labeled salient image coordinates (point features). - - self.features = features - - - def downsample(self, factor=2): - """ - Down samples the RegisterData by a user defined factor. The ndimage - zoom function is used to interpolate the image, with a scale defined - as 1/factor. - - Spacing is used to infer the scale difference between images - defining - the size of a pixel in arbitrary units (atm). - - Parameters - ---------- - factor: float (optional) - The scaling factor which is applied to image data and coordinates. - - Returns - ------- - scaled: nd-array - The parameter update vector. - """ - - resampled = nd.zoom(self.data, 1. / factor) - - # TODO: features need to be scaled also. - return RegisterData(resampled, spacing=factor) - - def smooth(self, variance): - """ - Smooth feature data in place. - - Parameters - ---------- - variance: float - Variance of the Gaussian kernel. - - See also - -------- - register.Register.smooth - """ - - self.data = _smooth(self.data, variance) - - -class optStep(): - """ - A container class for optimization steps. - - Attributes - ---------- - warpedImage: nd-array - Deformed image. - warp: nd-array - Estimated deformation field. - grid: nd-array - Grid coordinates in tensor form. - error: float - Normalised fitting error. - p: nd-array - Model parameters. - deltaP: nd-array - Model parameter update vector. - decreasing: boolean. - State of the error function at this point. - """ - - def __init__( - self, - warpedImage=None, - warp=None, - grid=None, - error=None, - p=None, - deltaP=None, - decreasing=None, - template=None, - image=None, - displacement=None - ): - - self.warpedImage = warpedImage - self.warp = warp - self.grid = grid - self.error = error - self.p = p - self.deltaP = deltaP - self.decreasing = decreasing - self.template = template - self.image = image - self.displacement = displacement - - -class Register(object): - """ - A registration class for estimating the deformation model parameters that - best solve: - - | :math:`f( W(I;p), T )` - | - | where: - | :math:`f` : is a similarity metric. - | :math:`W(x;p)`: is a deformation model (defined by the parameter set p). - | :math:`I` : is an input image (to be deformed). - | :math:`T` : is a template (which is a deformed version of the input). - - Notes: - ------ - - Solved using a modified gradient descent algorithm. - - .. [0] Levernberg-Marquardt algorithm, - http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm - - Attributes - ---------- - model: class - A `deformation` model class definition. - metric: class - A `similarity` metric class definition. - sampler: class - A `sampler` class definition. - """ - - - - MAX_ITER = 200 - MAX_BAD = 5 - - def __init__(self, model, metric, sampler): - - self.model = model - self.metric = metric - self.sampler = sampler - - def __deltaP(self, J, e, alpha, p=None): - """ - Computes the parameter update. - - Parameters - ---------- - J: nd-array - The (dE/dP) the relationship between image differences and model - parameters. - e: float - The evaluated similarity metric. - alpha: float - A dampening factor. - p: nd-array or list of floats, optional - - Returns - ------- - deltaP: nd-array - The parameter update vector. - """ - - H = np.dot(J.T, J) - - H += np.diag(alpha*np.diagonal(H)) - - return np.dot( np.linalg.inv(H), np.dot(J.T, e)) - - def __dampening(self, alpha, decreasing): - """ - Computes the adjusted dampening factor. - - Parameters - ---------- - alpha: float - The current dampening factor. - decreasing: boolean - Conditional on the decreasing error function. - - Returns - ------- - alpha: float - The adjusted dampening factor. - """ - return alpha / 10. if decreasing else alpha * 10. - - def register(self, - image, - template, - p=None, - alpha=None, - displacement=None, - plotCB=None, - verbose=False): - """ - Computes the registration between the image and template. - - Parameters - ---------- - image: nd-array - The floating image. - template: nd-array - The target image. - p: list (or nd-array), optional. - First guess at fitting parameters. - displacement: nd-array, optional. - A displacement field estimate. - alpha: float - The dampening factor. - plotCB: function, optional - A plotting function. - verbose: boolean - A debug flag for text status updates. - - Returns - ------- - p: nd-array. - Model parameters. - warp: nd-array. - (inverse) Warp field estimate. - warpedImage: nd-array - The re-sampled image. - error: float - Fitting error. - """ - - #TODO: Determine the common coordinate system. - if image.coords.spacing != template.coords.spacing: - raise ValueError('Coordinate systems differ.') - - # Initialize the models, metric and sampler. - model = self.model(image.coords) - sampler = self.sampler(image.coords) - metric = self.metric() - - if displacement is not None: - - # Account for difference warp resolutions. - scale = ( - (image.data.shape[0] * 1.) / displacement.shape[1], - (image.data.shape[1] * 1.) / displacement.shape[2], - ) - - # Scale the displacement field and estimate the model parameters, - # refer to test_CubicSpline_estimate - scaledDisplacement = np.array([ - nd.zoom(displacement[0], scale), - nd.zoom(displacement[1], scale) - ]) * scale[0] - - # Estimate p, using the displacement field. - p = model.estimate(-1.0*scaledDisplacement) - - p = model.identity if p is None else p - deltaP = np.zeros_like(p) - - # Dampening factor. - alpha = alpha if alpha is not None else 1e-4 - - # Variables used to implement a back-tracking algorithm. - search = [] - badSteps = 0 - bestStep = None - - for itteration in range(0,self.MAX_ITER): - - # Compute the inverse "warp" field. - warp = model.warp(p) - - # Sample the image using the inverse warp. - warpedImage = _smooth( - sampler.f(image.data, warp).reshape(image.data.shape), - 0.50, - ) - - # Evaluate the error metric. - e = metric.error(warpedImage, template.data) - - # Cache the optimization step. - searchStep = optStep( - error=np.abs(e).sum()/np.prod(image.data.shape), - p=p.copy(), - deltaP=deltaP.copy(), - grid=image.coords.tensor.copy(), - warp=warp.copy(), - displacement=model.transform(p), - warpedImage=warpedImage.copy(), - template=template.data, - image=image.data, - decreasing=True - ) - - # Update the current best step. - bestStep = searchStep if bestStep is None else bestStep - - if verbose: - print ('{0}\n' - 'iteration : {1} \n' - 'parameters : {2} \n' - 'error : {3} \n' - '{0}' - ).format( - '='*80, - itteration, - ' '.join( '{0:3.2f}'.format(param) for param in searchStep.p), - searchStep.error - ) - - # Append the search step to the search. - search.append(searchStep) - - if len(search) > 1: - - searchStep.decreasing = (searchStep.error < bestStep.error) - - alpha = self.__dampening( - alpha, - searchStep.decreasing - ) - - if searchStep.decreasing: - - bestStep = searchStep - - if plotCB is not None: - plotCB(image.data, - template.data, - warpedImage, - image.coords.tensor, - warp, - '{0}:{1}'.format(model.MODEL, itteration) - ) - else: - - badSteps += 1 - - if badSteps > self.MAX_BAD: - if verbose: - print ('Optimization break, maximum number ' - 'of bad iterations exceeded.') - break - - # Restore the parameters from the previous best iteration. - p = bestStep.p.copy() - - # Computes the derivative of the error with respect to model - # parameters. - - J = metric.jacobian(model, warpedImage, p) - - # Compute the parameter update vector. - deltaP = self.__deltaP( - J, - e, - alpha, - p - ) - - # Evaluate stopping condition: - if np.dot(deltaP.T, deltaP) < 1e-4: - break - - # Update the estimated parameters. - p += deltaP - - return bestStep, search - -class KybicRegister(Register): - """ - Variant of LM algorithm as described by: - - Kybic, J. and Unser, M. (2003). Fast parametric elastic image - registration. IEEE Transactions on Image Processing, 12(11), 1427-1442. - """ - - def __init__(self, model, metric, sampler): - Register.__init__(self, model, metric, sampler) - - def __deltaP(self, J, e, alpha, p): - """ - Computes the parameter update. - - Parameters - ---------- - J: nd-array - The (dE/dP) the relationship between image differences and model - parameters. - e: float - The evaluated similarity metric. - alpha: float - A dampening factor. - p: nd-array or list of floats, optional - - Returns - ------- - deltaP: nd-array - The parameter update vector. - """ - - H = np.dot(J.T, J) - - H += np.diag(alpha*np.diagonal(H)) - - return np.dot( np.linalg.inv(H), np.dot(J.T, e)) - alpha*p - - def __dampening(self, alpha, decreasing): - """ - Computes the adjusted dampening factor. - - Parameters - ---------- - alpha: float - The current dampening factor. - decreasing: boolean - Conditional on the decreasing error function. - - Returns - ------- - alpha: float - The adjusted dampening factor. - """ - return alpha - - -class FeatureRegister(): - """ - A registration class for estimating the deformation model parameters that - best solve: - - | :math:`\arg\min_{p} | f(p_0) - p_1 |` - | - | where: - | :math:`f` : is a transformation function. - | :math:`p_0` : image features. - | :math:`p_1` : template features. - - Notes: - ------ - - Solved using linear algebra - does not consider pixel intensities - - Attributes - ---------- - model: class - A `deformation` model class definition. - sampler: class - A `sampler` class definition. - """ - - def __init__(self, model, sampler): - - self.model = model - self.sampler = sampler - - def register(self, image, template): - """ - Computes the registration using only image (point) features. - - Parameters - ---------- - image: RegisterData - The floating registration data. - template: RegisterData - The target registration data. - model: Model - The deformation model. - - Returns - ------- - p: nd-array. - Model parameters. - warp: nd-array. - Warp field estimate. - warpedImage: nd-array - The re-sampled image. - error: float - Fitting error. - """ - - # Initialize the models, metric and sampler. - model = self.model(image.coords) - sampler = self.sampler(image.coords) - - # Form corresponding point sets. - imagePoints = [] - templatePoints = [] - - for id, point in image.features['points'].items(): - if id in template.features['points']: - imagePoints.append(point) - templatePoints.append(template.features['points'][id]) - #print '{} -> {}'.format(imagePoints[-1], templatePoints[-1]) - - if not imagePoints or not templatePoints: - raise ValueError('Requires corresponding features to register.') - - # Note the inverse warp is estimated here. - p, error = model.fit( - np.array(templatePoints), - np.array(imagePoints) - ) - - warp = model.warp(p) - - # Sample the image using the inverse warp. - warpedImage = sampler.f(image.data, warp).reshape(image.data.shape) - - return p, warp, warpedImage, error diff --git a/register/samplers/__init__.py b/register/samplers/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/samplers/build/lib.macosx-10.6-x86_64-2.7/libsampler.so b/register/samplers/build/lib.macosx-10.6-x86_64-2.7/libsampler.so deleted file mode 100755 index 4f248f0deb0df907b62a956024091b4efaad1008..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 9528 zcmeHNTWnNC7@lnxSVHBrL=+>JRf{1J_JU2>ST$YP)w5xX)lwuVj;Fg@x^Zvpp4FCM z1i8fZu(0t(jXs!oi7}0k*ccy-wh)zSB88~oL5Z43Y@&%mBtq2n`_GxPoNfz=Jo#|` zWX{YtGygaL{4=weotYDte*I?-V~$3~mXrNQ$Nj9 zArly!Y?cUAgL5@8s8Y(1dOly1k>A;1#DrnSWpT1ymjUq=hSGi66fda7KL2$4a`TPa z3xpD3bfYFN!8cpvfDc4`Ad$)?^2n>V?^(g;5%o=6Bb;VmG4oV16X(fHe|84GBZ6;0 z5D=cy-x{%z;5;$F`}0aVVdAOdu}CuVe-?U#X-rL-$Jpa$nqnH)=k#$61vdsQCFXcu zA*OL9Uei|2(^=|dQsw={ObmXTnUCg{F5tD=Mf1%noM$km#?a=}{4_q7uxpV}A}pqi zox%5T#e6N$4A)8NBuSD2sOZYU!gz=A@%(Y2tMmvGwV`ZpHJf;BXAy#Y*-%- zMa(tkMka~iV{RR>04*vA*=6Or8JK*sX}9;~lTBwX{B+(sHu9eD<69ZqlkU0fSa++h%1OmpZZSm+rYOWrcVHwm_MwWFl31c)c#TgVG>f)Vxk|l9gp)&`%l;n z+Z*Ju?;9hr30XRP6jpx_(Yjso*pH2JY#do9<;E~NNTR$x?YNEb0eB0kNrUHeIb`Ffn|aSBjwOSb}l5Xbh;v%vlZ|YCOan1NzG1I zs3SZq9cYKwWvx*+NIHDXu(b=eUVrhN)Y8;YYk%jVYIsDBg-2lUXrzNZ4aOTLzr>*X zM;{438r&LuEXc=?T^>LhKMebKaOL>qLf!UV_@WDj|EkO1K~Z)2u>O4)-)MQsAwhfn0R|f7y^~tlaa9|mwzNHjVOrAmtHP~VVOeofhKzN0~h)7z@{|KD`6dYm?O= z1GWcj57-{CJz#sl_JHjH+XJ=-{$CzwxrOBx_2AH?lXJdpaFFBryO7N&siZo@pY*W@ zx4S9}I&lFPH!jHgaHow6^2&BMX>8fhqXv@M0`F4_iM9xCf!F!NPqGTzzGCL$kSO*o zp!Q;e3RmzhAgMw{h=XSmNwzQ;mk@GI(H>L6I1Y`-+)(DuL%83e9dz>^48W z2vO8zOn=;v#DdrXl+ADjk2>O*pZQ`;$G&6|3in4JSR0IT9;r=V%x^W;&wb2IFH79U z%aQ_d_o&eT=gBO+nqZdtS6EFmEi-)`OF`0^caYvu?q?K5LKB1&%9xn4| VClgold1c5mfNEj{>&Fc9=3mdCysiKM diff --git a/register/samplers/build/temp.macosx-10.6-x86_64-2.7/libsampler.o b/register/samplers/build/temp.macosx-10.6-x86_64-2.7/libsampler.o deleted file mode 100644 index 0f2d69ccfb889956229d01eb92e5e21546d776e6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 48620 zcmeIb33!#&wKu-^IY~GKLIQ+{BuGd=$p9n}LKq7|G*LhVQHsd% zDo(Xpv1--Y!M0j$uTxtGw5>y1hql&QJKWB#z4g{>Cu;}4-*2t`zVA6F!K?SV&;R*< z&-3Mhv)^~Ewbx#I?X}mQ_xt=?ul;M3F*a=QgZ~U~E5J`B{^RdR{BPm=VgBqF6miI- zaqr-1_*iB}ckk&X3V}cJ%VYSUF(Cr-%zg&mUo_g8=ospb zb`8dQ697h|t4~|G;W!U%UUqz*c94S+&D#s$N51$Iw))8U$}id1)lXP1zuO%3bvrc1 znWHs)9sFD#N8}gl=on1w9L}$)Pvf5ex%l^|4BaSXWcgx)vA*s^YM6X`kp}+5WPWc|c+aTXNl}As1 zT$0R{FY*Ban_?O){+>Qx!$aV=LLW!W@4!$;A1o%7%kPN28tEtv%=a&H4?%D^e(5~} z!}%?6{BCub`}4cxm+YekJBR7-#S2uB4NjixF4KoYUw&u(rP6~!n490DFVy(k>NWhe zFm_6-aXBn|!av`>o^jvdeD^pWXgd(Pto5$;%U)`~@X1Kqfen$% zHr?BH@uA|YHX;{Kwcf`>)?Rt;jL6z6SI;OxjN96-TonP=gI&d^J=1pacZ#nX56U(1 zGwy3YK#Bt=Mp_%c)A;tQHy~zL@idjMl`J0zH!_lgt~di|g|u?FaH^;PvaK&A#wPaYe4unOZ z;(NZTvaW@!rGI;)c>LIqpZ>X5zukIYTfFr?2z;o0i8&7(%U*r%UvIzNcwu+t%3KC1XcLu7NK3H$wyC)#t_TJ&*WtC6|Qq=p~8@JUZj*75VY+ zuyl%+x;7GBlfzfK@Bt@JF26iFrVB@1QE+AcRe4v3%{AyRs3VQojll=ocO?inY$^K z7);IW>xy-ECu4ncS0_?C()|N-2lk}5_xH^nOmrtx>BQh%|G`C0Nx#mC4Pk6SD%p)D zrZT=gHi*~AP+t}xl%Lw(KbWpe_N7fIG)iB58iD5JGYChcLkQ>3i>6K9?AMS&C~s62 zAt*r5+(d#oue(xygjHlz0LY5#k&-HJMzq1~&s#GD7WvU=Z$}EgKG9{ufz{86Ax}f@nI{Z6@cNi-0Wb9*FjL zb)^z%ctKHRwJGcx2Fok#-gnGVjSY6=ei3Q_XWRBUG=8(w0@L;~S< zQ47pT*Nae803RIclML>-1k3_3w2~%uUp93d-Lf5or>KvBJe}+nE;|jCtL{HB1AioiGf%=VG0*1Vq;-nc7(+WKLzmY2ul>P ziJah}kl$3FDO{@1Q%N2`mnrl#kvd1B<%-)ZxSYh_r?~gKc-aE1P~;XDF&nv3k*9kx z16ht$%vLX2fLW!OGfYv0PvcVQPE^-#pg&hUX32M8n#E(5>f5q62pqfchzSf@rZ6bl zJXL}%izPT?s|3+K5_DY8U|eV3BnF8R3A!gruzj`!$yN!@IzxgTLlX2{DM9ZiCFpxv zg8u)IAXQXGX6XeI?A$8Bu6+{hzEgs89+%()uSjt2s}h_SI)d>oI6{I88zi`BodlP} zCHUZZ5?uaC2@c#R!4=O+aQ}Zu@IW#8i^by}Y?0uhQziJ)UIxWWA7yS!j(ZuEt+?d) zKT6OxZYtl}Ya}>fwFD=gDZ!d^C0KhWD=yhEgTeS?S4gmG zl?1B~G1>9Q9Vc(cZY&Zt9lc=8Gk`Qcf z$5%^mM@)iG?2+J;H%suTFH7+0?@RERze(`9@+!u=`&bF?>6G9L=S%Q!ACus|LlWHo zV+kI3OM(ZFfSW6x@KBQkUpiTWFP|&HSFV%b;Rhsmp-)M$>#Gv%{=Nj~{6T{A#=`X%PuyE0!M+nEIKNkd3$Bsi!jDOC(N`q6_<0E~ zd0m1JPN-w_OPeIPe7yt*1|@jzcOuyf$6@4NJn`FCF>e#UyG7n!{H(lv?`om^{#_Eh z91+|Pz9w(4{91w^)&Mx(&3c&`jT8F%(JG@~kv4WhtOxCT3@6%eb|ShQu{2tu zzW%}9Xezy@Ct>m?-8E$lI_eIX?CVUR``Mf9L)Xmg^T4r=R06G*Plkp=0M<`mt3Ae#y_x7~eFRW7WBy1$vm5iJH!V+VSUC?U+eSXH`t9E|&3wM1_nf<~& zwIeyO3onwu5T<`RhPuptO=W0cATbyhPMTD#^Q@s%x+|IJ=`{CB;!H>1zafH&r7Pef z5he+JbUA%Q6Mda#|NcCqYC35k(T8a))P=Gg%g~vkurUnPgfXjuFz873sd^BK{QX1d z^pY0SH<-Djz_bR`#;J8;&hhT{_*8U)c{)?LP*FQ^>`Z zGio^_l=IzG*?#(vf48fq0Baxe2~D*<0es4SVKjp-Ja z6pgZpQ<;s8DB!_FSG<3yFWqXGHygGy0Zi?|=&3WhE1BL7AKIDhH)Z7`ahTUI5^<}! zO<7I1l%lLGWbFgU&Xj9Kwj=Eo5kQXsR3TGD1yDor5gS@gWpK(82{tX405fa^Zr&=v zmR%BTy;g!V?v`NNvl2vqE5VuJ(-z=tmzt5KGhln7qkP7ERW=jHlbtr^6V{m}6X9XKLF&fbQ3c zF_!Ku%vE53#mcibx@BW}O^Xj;Ia^a>!!YwBWxQnVZ;*}}5}(k)e+@w`d--+tkr2j9 zCmM8Ftu2bMpC2Ab4|XPH$U5a9Eg}-h$D~A!{Z$#*oT2u`ZhpD&6_z!2Ey|$0?*YWO zmYsiXxP*d0Nc3CulFED8o>|LrK|o>ZZ8g^J%FjGv97f_oz!8IpjEL>18|5d{UaftQ zaQTsnA*_ocj{rJ5?9F9V9QP`u4tHTph8{x2X&kU)CL={(TT;<621{UtS^^buzR38U zgJ@IHNv@clfR1j%6Js6MzNXb2thgvv%4h?o;$lC}psBdTMJSB!G8G@xm$8LWZICK1 z{erM70U(K9NX2FPJ^{&hnTpGQ2(d%elMACg9bKCK0p&2O5HsfeJCcZUmDbY+KzjPo z%2ZrE06{|4ErqCRu})KQ&DHY0q7Wn0#12#Oq0h^E8zdS;SHM(U`>X)g7efMX_Pq9>&99GJQ(udtMN#U0~e7%=>~ zb~$<)6@Rn#7|@jS2b@V?WB0=ySG>iEJ zv+4luS`L^4wJMiJ z%P~(yp8X~|=*$&rKdX08r5o%$7#&oMwpyn-&Q=xkyU{3-|5mV1pNGTz^ zBvMQvlZ8eg9SW%+B$?{OH2c)Jf+`8>?oDH6*i7xzZr3cr&O!IUOij=tLe&jyi?gPm zkJ{FW(#9llcQ$<5TxL4#X22_tyPw*a5p6V+D^H$^?jq)=nfb~M9|J<#yJ&B0x2fEy z?|H<>22ABC`X0f%>=ZDSn;0EERuuIQp~-M`118^F`_hAZOy#N1AyFhV&y-YdPHHlJ zZqdE+{hCuI?SZzCMqAfvJSndY+#_NtPhS8Cc+EGZm1o@RcyYT-%22JITdjJeqKg`qsmGPgc>~YYN37Y%~A_jze!9Vrv>8FB<9X5&r{4lh{;0i zeN3{k|41#5^>~x=m}cfs6-m5nJ?8R6IbUH*^e{yc(nXt4?tdx$LU~)n}+6mdYP?6pG0w-BX+odGH(^O4uRr0 zaqMu%GY3`nrm-H5v&jzk6^!Gi*Vqa)Z5oeL%ntV@bmpfouzx}IRvhP@9WH@ZW%?5P z7Zjl4I2o-tj1tggT4*0YBQf>d$KdtptF5iE@jO8G4uR5I{+53YZcKi#K$+LsKSZ;v zw-+JG)C*p~FY(siL{k310yxWEV}Do$ii;KdC1UvlEFm>^4FsP00UbenpIH6?`(p&P z%77*}^*kkcizNI3i9riit6a`Ozh~-xCTlI{$d#nVz735{JzwFRG*@`N{T7?$Dt&vh7KVs)^h$`bz5VvLrGQb zR5&N(wL?*3SE3tQwd*B7tmSXH{IG14qA2nsKtp-?-g0V<9YRg0+O4eeDnUpJkKbY# zeAON;+v!BgkMtUfBF6wwOX?c?3N%@Dj?%T0PWlFX2eKZzs``K;yM5$zJk1C8c9de( zxr)7b7WKsA5$Ys7Mp0gN#35G8sXwa{4Z{c9C*gO%mZG{w+wh#$JLtR&}w8$<>W4 z<&K62sJcXvT=Nhlq>Q18R$Z#ea}h*&3)|@sno64qaMfi>!UYm77Y}`2bu`)iqiGTN-UpQJsSv&AeZ)Y`Kn<#g=`os*k8tTwTgw zTW|Nx1oGetcr+#8x>K&^#5dH~E{qncZd5X^O9>xH_BM0SpD^<*ng7J_H?A-KMEeha{?7c_Ogf)hVj7#*Ur>;@g#&>tjKJ zK%PlE{iOQ(n2+QWsuM-1vE3y3(nSDS%U_KzB%!=KnbihWAE`dALUX|_LuesIag*Xs zcAu@~T3o=#bDCdAW~)BSjUQG*9eDG?4|A9qZC{|R0=VAC*LV_*!QX& z&v|mtW1L2x@C8M1tx*IdVs0J~*D6o0KC0y831phpEsmzM` zIv7-aT#;PA%!>IzHdCKduiJzZ`c?Is2vX%>m<+5q2 z-u9v}C*Ic%^>K;O+HZO|B?wjvn8lak?z30pENdXamd4H7|ykr$jlfQs|(K_>Q9@QZb)pf|Et%`bfe)s z6Yfmz#9Zyn&h;8?ktvzk#l>wGty_o4F4}f|ih%O!r9utS7z!&A+?T|2#8^Y_GAjp{S=4(i!jZUWHBIEsx(VW#p~ z%^SA{NF>j*{RK8Q^K zy)Ob1=h|xQ18^ZTzo6RVPJ|$mnQNbB34YO9?oB9v4LN`Vd4*Z{w-UHMEjn&k&}_)y zvc~>_1ovw>aT`N`^r}J}?sw)R+JJC}LoP+V-CPKYCz`-ob;S)3LIov^N5en!tEwjM zh;XfHC7zt~-ozY!O=EH=M9`-}^qbMOnE7dSaNHjuzD*rmjV;Ot@z1Wl?i!jDfqfL`VS0MCwXY+MDFj&$ zc$~F`vpryWDNw+iJqqL}T07InOp!u@Rqewt{Frrrg%{Ntbv#L?kz+YC;8C3X-7&_n zMp&B+C_btbRpY?p)Kkz$S02~+$&}O9q`I5s_I8q;d94gZ|1BV49HB#6PnUwE^LZ4s zn%N6}T8RI5X8NpKQCR5yg_&dhRP-VyNDuTpWL*-Of}GXzu}B`5 z-A5hIYCIp{BM63fwR7#W)$kde)}93L6ym5}0d^t(K2W!^tgxR!|GTCKDGI4mEmGjd zm?e=b@mtHsTCV)oK8X}DoP#=-aTiZyR?EU;jak`;1wK<+VlP4_39$B9EW3Bc(wJcz z0;be53d9Rfh;HC^PwuTitBhx3SJ^69cj*!KZm6633|l*9E@WRyPR60j?H^zV4n3FX zL}RQsXbnbr33(I8xEtKfz}^i~d$Sej6- zvJa#DN>|&J&W=Eh>1U@6Wl3tJJ&`Qe*vnvcEF2IX9}v1^*MhWklRXMF!0!NF-5Swq zx2x?*XopI-*&SdCnyW{oc?R{PwA(%i1%l?^d>Rp)>(|KQCX{~ZIrbOI=leblRa(v; zk7lo(g^@Y<{N1NHi9IF$cr+K-$Kj7bGp*5+S?U^pJeo`F{SX~At9%;Cms8{?z@YSU zyBA#PFhJKuIuTh^xZ0jI7Ax;xwL76}psxYinUzoEDd{xKq}APTFCHh9pYthAtG%2y z@KO6GG_xDk_NLn?VdEp#elMV$E6}>fto?%z)LQ5;4_Z5Vo-e}t!3#B)oa(-A?d)t! zMm8MH3)Y@`IGU_vu({=&fY7&uESjHNdt(lbY<+WSy=v_v0opeLWij2(nuXl|9xH1< z@QIvbX|%iGLQ4N+pNAr&$f33W@M+2=%Z}t9_D4R+wE4bPq--6@sL&=f0LtaqfJDr( z(O$$lKQVL<%9h4w?U_D_aoD3m9LUh@tldkjx|~M4uM%;O4(%W1#l1N*u3#HO$G~f7 z+%IRw{VA$=>Eh5XnD{K1m$ffv#uesk62IUWoSN1&j=-vuC_nxGqE=l8xMUBfnxhIv240Zw}&1<@`}A8z*gC< zs7Q5vp?aAJy_XpEy?(i4M-*vFFAOb(gFw$On+BqPl>YfZXj`$Dz-yTatSdnRzqCsw zfw2p+5&-WuWVrNGp?`xtSbMaO2HV*afw?=xEdq3b)~+Vj)5sr)eK_<9B+4e%+HN1K z)wO3T<#vmMM(Y}S_ETKE-wt@G6%TQyzmw1 zDyqT$F(7efFC<0bUqh2hQr7Hi!kN90l!Xt##FS)FK;q2)LYQvd3HB54MoO|ZAaMx^ z$?5hfUzW?VNiuv(qZMLoz6g2zz7P=mz5_);y3%9Am%^|$-50az97)|#_9b7me*`3% zwzRJltPJIMnw0MTY&v96NbU)b^;!KmAj!<&?Wo_S-wJtI=`#z_bOo)2SPtRLeoI6v~H=N^}sHItgo_XL!P>Wp}#qYh$!1eiV_sY zarRNDsqRw;QT3=1#JCtS)I4bo+K+3r_psrAJoI-o-r6nrA{k5bs>29hAhgm~!!GOo zEQ`um0|+EwX+ss~qAyhXkNov0Qc8kOu+IzFt?c*|6jV+2utm$flz=Tj-6}0AxCik# z7KKhSvG<}1)jc120C82y%|2EqLjoxqEs0BO3XX&OQ{w#rF>PCvfFy-(VZrZ_veMld z&;^=Qy1HM7-b9n7bT0&SqDm$Xe#=m;N;eg3Ls==^n*m*Bx}>XnGxQWHmeS2!?kQq; zW`yb41uM`)QR1~D5WB?dio%aV6r~#q=rV=83aOS}UGUeVJi1!~x=bNSS9etSLq5N6 z1a!V~JWUGI&lQyTIeKFRVhAa8j}{c5?N%YDywA^hrjQHZl}o=>u-Mo5u>oDCkff`d z9X=hM0Oi*o&A^Nj#{I1j|Fs@LcYRIOV9SGvxCE>p;V!5B+dk7_}c zP`WDux=bNSSNHz#GQVWL6wnP5Qkcd@9rV@yvk{0Pq%iFZzjCCP_~?~>&NGE91VQQf zqaM#Q=HpWL8w0w`iv2a3#fL}5;j&&d^bx!l#`!8{21$~n?t<`*a3e}`en8^0@|txa zzBYV)xkvoa;fS-Eb%}OI_>Spbv>yi1s4Gv5FNRyJM-qW}PSQXOA+A=LzZ%{O-K)$a z5^7yDwJgS>$g$v8OJx;B^!HbIUJ}t?B9{n!*$)7oTzcJS0V%|M9k!ieGovZJo?N(6 z16VwUW7)~{9#bmo>e!GsmIr1^Z_x8HC1x49PAwuXMAECv7B~j$I^yQ)|Z%+9Yp8}Fp zjso&W55t4d-C&Xg5^93d?0RkQcQBBu8y{-U(-}=1H33~rjA`%VUfTNYa#E?8heUf4O@y_YJL@T#IpIY}Sf9M9 z&^U-av2Un1F^KJ<^=EBDtx>wr36~!Qx*b14heRQHCvYEHeb0ZOpcE20;jR!Mz0ZMC zAtfi=kJ?k;w{4sUscgq8w5cD+6lgv`o;>fXI0=qhCb5&aO2yPqwaf=6uz83|c4AVr zN}I`{k!PkHuGKJtJ9O3E)8BzpjrG-8x6+L+fE1rWh$07hvW!8THPVhJw*DA8YRH6VcvbN6$mbU}tVXn6J(f}j= z$HlqD+PeY_dD%8X=_D0?jp$uKvq70^!s6t2mqytp9p5ss$t9_`~s@ z%qm`35Br)~OXTLZm07+09%^&_ z>tFY$jLyxWXtXcVq$b+uz~l+oGp0g8qS4pdM-Z|8NV@{PO%->Y>UAl~z&=f0qV95* zzd6ZXPt$v!-Hh4{{7H7wEW96UzdB0ZSJ^Ycr@qxrK(SF%peaB^U8`Gff5F7sYzQtF z7{^)E8t4|aZhtFkT?^tBMV8~YmStBIq3lZ{ClSEMB6&P{J_?|u{xrJ<@Yo5)Y>4i} zR)HB^xP2pSxRZbmd?THF{Z@A>G7`o<_jpt*u9@A!Xi^V6P<=GZHWvspIfUrNYBkp& zJ*a5S$5}f&>`Bn^F7N?}tklF4Bd{4}=uU@8dxh$7ou%d1@3c5}*wzmHpG&s7s(XrE zRtX#`0<#N(vwQ+)eyi;j@L~1)?2+(9Ah-$y4ElF2I}zO6#w;M->QaMu2kpzKwRbf?&<}hw*4O$uxWY9rpGa#@q-9ttSHyW1zLk zwtbd@_-|`!yYXBUp{|xhC~#5aDnM!#z{e$mJarL_C{q7Ld+Ld3etDp_{s}pqL=P$F zQtQ8-Rh09Rog<$~U9EyXZEry1&8&mr6mqrYTk&hwSvs`(Z!nWHb8yoYhjkCT#7CKG z`xLZVj|0IqUe9t$HyJWFw@MZ|@LJ2qq6i;PW%Zd={{uUF9daX<|9fZoCd`FQG3S)2 ze@!UR2g-ordKny7tp87!{SraNWAW|SsbBvG7h$qM|c7uo;C z{-iJf%^m~LU+pqrmTQX6_};SL1~X{+k~Js+X}mn+`n=1Pup+akH_-= z9ETryk~}Vn{1Sko$POXn(fXnY&ojTD;Frt3aS-rV&@Xe1hj;kgZ$GyD9mZN_><*MB z=}<*a=hY8oVS7KiD+&sou9rZR71%PQpuoJ-c^yRAXgj{h14Xtppi-2Lvr(u?sY*_l z+aby(*ndY^D6H~y-V;$a$u{^gSDdaFMwA_CFYuGyzjZhGlvUV+m}=3e#%_HI^`@-S z-U%f;SZFIZXqU~fUqQwlEbk2O7KyTBY(1*2gGbJA_eqr1+S|~-aNv?NyqhJ;X4`Kb z&L9v@j{KmOn!;WeaUJFbXmCZa*xu(0pJz{Z zcR7@;u~jhde|HLZw>p%aWD(s#p>ub?L)j^Ip%2bG-Mbs2>@;h}demkc@1-a^OMSV2 z$H06*Tg@*unQgpeq6{kuS=SKE2h<}8ZNps_F&tAXOXW~CE5akp+wxa9`?8((e0VHn zAKCH^g7(_CZ4gkoJT| zcL$I0Xlu&~L!U=|LmF|P2=g#4NTqB>=*daAXBMSk?fJwOf}QoSb3!-492ENrKW;ez zAo&$}FAv=VYgLk`MkYBabaa(R@`r2^WmRPtq288l3w;|N$J%kH_?mDjJv|fXvqHZ^ zf;276c96sqRj`AhA3-OIJt<&707U&LoKo5T(2rs0O46IjK(LpFehNpY*lUI{I2ZM= z?50ow`nWo?1&u0A03*R5J_jlg;7@E2z8Km7ho>a}nJK`hN`ZYi z^f?r`4ywk&SX}7sZuls)wG?PcMBCjG@U2wb=jQmMDt%#W>VO%)v!J^7G9*WE|>Y6lC?~U-Dz|UKvHz*lVdS zcZTc<4_8~Ae!m(lFC`0WOEv?} z9y7N29)Ph_f06pcOrQ_1L$ zGQ)mCyA!>K$uAfE1bIsfs9dJp{U2By;W)d28_i-J$(@ZR-<-;o%}B#}p}{TMvF8|@WPV8bn;KY_o_Q68O%{T)iI;r387iu)=?u&asHL>sMO?+l#)7T1d^x_zv6 zj>Ks3SL7Yd=OD1$qFtjB;Gh#I~WT3zW$|4ovJ^ySxS2Dxdc z;gQfbl%S^Zp8*#~vWz--JhXGNm(X}b&E~S@7?6BDR5Hcm(nJ!;s+Z95dBVnW^%an4G@pzvTk@1A7LD0Ka4Dd%iff?$k>}v;Nc_mjh7mb zWl@AKZuC0btQ@ebpGStxKRMl)4q@2h7=}B9<08i~+#xI%JC@-%p!6`?;Q9m1^IL(bw*_eE!pY4?;6w_~FJEBfu{CCJ-JfQg+9y9WHUVzEU-<_yc zn!j5Eimdz=R?WMn z?&X~dSrT~?(Q5fv6rlwyiF}^`K6bW3Io3|xW=!JtQkezr%Lf=)(QkJHcu@uh2-Q zlek0!7nfsizMOhm0}@~c29y1`$7+v*_A_GgukdZwzb*e>Wp*tBd)5jzEDF$%QSkXE zMEwp)?U2Jmn@lI|5?4Y+oi_Mjd!+iC@U{tfGqtp>S{6?yE2*)IJz`OvPfW0*W0ghA zmqWKZnL@SYqLAsNNu$T_QE0MB^r%>jy{3~@8Z%xVUQ^1z(=4n#45 zSTnW7v@%|hI@E((SWEbJc4rdDE139BBzAV^S;kJ}+hcee#C5H>+D)5_-7{x7#R6G%Z0WK_x$MJGhofXvyt*-@>+>r>X1N?x^$@b z0J4)BU*~Aw??PMg)|prh_z2nWrpi4;r%$l!0hjik=;}mQ3|ZyvH*4J8A{TwVf~9c| zE~ZWON4Mip`19Exe9#b`Muy?*#;^BL{z zp=2T*T^Cyy#Jjf?@&0)DcnyrFH|w^V8(d-{6E-5A49YZK{O z-EN7!ucPIyTk+pT;VubTRlx;PvJBzJ*WDKHIzDc_?)Kb%r9KBe6!t&nMK8crz~pE7U1-R)=>45MvUluwMl z+CKdTM`QCa;st{7X;#L(+b9?ow+ZL>+E_= z_fm^?!CBD8&_>lKB~aQ zz~vD`y>ez1#h5|T>ZVhIc(~Yaa@7dYcpbBHPWE}R-pxeK`SxxMuqQW+bCSrF#gf1F z_?nJJG#e-cxNYrMyfbyfyQkjuo_K@$Fud<9R=F1ja^p)cTn4B3|D#BU$@4p;^gj}B z1PvHF-p-6n921U|M}D?%>;_wcq6y6~W5z5841MV`6OJq%vt&#=W`y#Nu;IvEygfh)> zixAmMu{dJ*I2rH~dCceI_(hH(uH&7E$EE0JOw7wlvvd~V6_~1;L|K-OrJ!Ya94EQS z$0pATNZnF650XrF^5qkNvE>~5zOf=&{?f5f5ys{Dpl!kfwrNtt3M4v?iB2U)Wqy>O zem>Kh7No`3N|#xe#PpmbNKzHRKDYu*kut&6S>X1GPNW%4s=fRVu86cxeBao4V<^u| zK&iWI=+RVxm+xa7$rvQUh>SKXfPonyHGZln6JM6vqb2)wbBt4vhcu=>h#-vpII}Ub zHghI9TbUb>`4kO44BGLbjR7>9;XFTej~v@VFyUyLB%Ar3=`5Iovf)d1^kxT!ijaRl zxrIL040>EcL#qr!DVRpVB(qt(*tf?e__UZ9?9z$m_=#FcuwQhd^W3`DD++pIXo9OuVSHDts&9>N6r zRKsb@G8~1q4}(=e$`gjgCiKLB)@6{;HI7zg^K5%wJmxGla#Z5ndgJ-W>v) zPy*si@X3Cm<)cQU*fvnf^aiZ-8#C)JHGWFKWzCqW6El+C6kx=>{De;pQXv1+ro)yt zPr~dOCJB5LcoDEAC;LmsQYoi9Re4HL=+@!1@NH)drv-Fd02P;^9HIe-4Ct9bMwmt4 zqGAEHZ;>5NKd21qPvd^%th&@W2UVkTjI+%|KFmqIY8M|_Q(4WrSwEq?rm^Gg1iPZZ zk#7f$YkERxuyFrJIFrPfvluTY(hj=o5L|js5F4{%(5C4NWNe?9l~R8|%cuiRYa-m) z!3(7r%yaLFn=4-tHCf=~{FdfyNNP_i3y}2^u$hbR#5h8Hb#!ibJU%zs+|oQ7 z+k<)%v-RtkjSZTlC&b}cjTCWrOLMfjDI=9(8JVm16=h||&*fbCbbQYDP@xxo(1sg$ zsGvSwX_~X0iq)I8@r=S*%8#q2XD9oXXJMdYp9)&X_AXo=hZuNHC3+I^G@ePLdtwy9 ztI>G{GLnmfx--$`9{C9?4~j3bI6OYlDUUL!lqY^jQXb^qSdR`F<0D3${`atu9&Ayk zD~a>=}mXCLw^ zb$GTtAukl37FmjpZ!B^bpDXd$an}i<5>(y15QrItf8lQ#P}}fh`yzKeofi)Oy9w1> zynlAzQoOxn$m4VS?ppq*@Yu-aaQcK$%@J3viA)Ob4z=0HRUrA9c`5vYxeH|9-Fw&F zk;Knppe-W@t(gHz)HMlZ^bBNc4GM}`_JR`q>E?GrfdRi zcW|g*`Xs01l6Y%L*7-b25^vSc&(HE`ZPZWjZX^bm#3p?bBDYS`))LE z+4nyu&CQn}<>x%IKl=P?CYASzC4P9>kLukk{;^CT9WeOap9yoHr{JeL{K(!HJbsx8 zbl))M2ayl+9cF&h$>*WpSMJON<$O1TyVQL0VEfuFIbizFw46$I>fnX&P8fa~OTVVf zy|=+PviP-EPc8g7jr;5rzl?CLGpzm-i+_%{Scn9q{&Zs{+aRck_rB)d`pWBDdFdjdGJN70$vFgY9EG~*L++~G9BEy4f07ZG!v ze}vsqT+Pqf@`NoxT*%LXVj~oHkmNwIg@%8~D%?|qyOOv+DfEFegLV9C@$44R4L2om z*P8tSfL=qu-7Erpi?QHn_$OapG?>FZabv1+9u~1Xqxc3%8lMwuU4LS{mzaz%4>s$0 zgVr4ILU-Vx)pX%IE^)=WL%<8&3cJQhEyQm5UVN&)b$#6DpIRVagKX_(a%)*e zk`z9vg*ZLlgOiJ`>lXQ!R%1;f_d@CR_3@3V*4B;lqU$y`MBBk!UKd1_3)Wn4a2Fqz zl+UT8koSkP>0uvz~jkpi|8l3-auQSDp}YfqcDJ1`ye<{^_^@l-rD(C3CY`0}o(ZDXoA z*_TQTrcoN%Ht5Qr8H_=Zh6~e_Ao9A6Q5;;<2R+jYLnT=!jwoTQk?!y5bY;Tmutnr5 zmUgHfhjJwlXdw7lXR~U1m@XFi>TgAJG=L0BPBCc3cUZk6wPc@R_4Biwa(ldWk$iP{ za44RJL2pnCly8E<2r~=_o*fmM%rXjNkt}1zv7p&Xi`ck5hKt0r5|9XiQfcAT{Psab z$Q2Ge;T{*Ep>>#A3{pu({HlhgsjEwdx9Y&f*cJ^8rM6@A6yMS6tjO(-*sJ3)mDO3Y zE4~)!A7?J(JEAFAgP4S01)5Vs=ztBs3tS>yEkF-8I;Y zaz*tNAI8gWMNrA4tO@$}R3K5RfuwQR>0n}T&?Hm3Tr6?kb4`L~H#Rv(HL!+4rkrha zNQ%vRI}_Oa(cPz7@$*DSz+Feqw)~37_)tVLrAw`TIgg;!d5vn+bg;MxB%SES{=Wp= zS%#z1*9-@6xPHOh0+k6_%4zh5XCBEFp`Ez7QQuQ&yHdn zwM9U)r3==>GYO5-WCmGv_Kr}+=F@$u5A(P>Ww)SuWo*}!+aO3N!pORRwlG+LA1MN4 z%k=b+P#Tc+QcK>AE;?NiXZVkW?wDk&*||&9(nY5z^$5Mz{@xCQQ56qPqfW5Z@|sGq zZh?7uZK8Ka{`a=?BfG4;a^gj%dOD#xx)^O|QD@=V(7WpdnQps~8aKRkNAn#;4cOJu z6YJX{W{^V7q6=dLFJ(PorbA2c0aHt468)7(XPP+KlaT$b2cmmn`e`j!>@(9=Ux<% zR-K}|+J|L}(d=Nqc0*E4sG!4)OZz05?O~k9*|VBKuQR55NRry(4Z-O_zCp-_lYQP( zo&Pi&h0a_VOfRSL(FswuuZHH`GI~H!doh>gP|V}z33ehq2ktciUfEyN!vCPA;hXlnLnFges<7Hv1dTo} zHbOBZQ}Et&tM2F~tN7YztP>wE^c<(@!mtz z;RNwc1Fy(!MCg5><3A)vTj2@^ItF(b=_E=Sxt=iVVX)B=NuuSD)=N@c=*C3tookZw zTU%SaP74cWc%42>4a-P@#Tj(TUE8Fre2blrlj^Oe7H4t?wGWeR(DQK|y;Vlqj2Eb6 zK8EA1yTvaAz{Bf_O@0qlPX7a^rLNk)mJ90Yi?b|b9O8?{3_=Ak1GTGR=9Z-uxL-H& zps7g}6@IBEx>WQY(P02VM~Z#HxaJrgoBly>z%M3%@{nC(Nu!6G-O$W*NRLaZ9pj_5 zOT}Y-uE*)>UrSKsm_<6(A$*L~tLYehB17!=M^SU>vAWU;c$-c&bsmXE&yJtXWD#J5 z!LApm`V!r-><&^&TV&UP!;9zjgRZRKL(<)FNZ1<}h#yqV`I=rJcNoAgd8RB~Yd1_r z>z9#NKhAs{oD<)Irmw{-N*0ukv=*Q#FtU-}wsg+frPZM;C5A?%RcY06Ej6w^ z;S9&-b+*h-zmCJLBnILAqbup2qT7?5`1C63PP%^xL#2&qnPk-{IuJ|4we^YHgAdmt zWhaj1FIrmD0+SlK?qh$L>LRNxSft$%xEiQkDz;>cIvJ(`8Mb9q%;7H5HKC{%IkwQO z_b{}=cG>t~J}DjV!Yr_uqbptMmj|tGP@Wkc;;_zueOG6&nZ2egQ)rxzdJ}wQiwc7U zd*f8E#TUQ9A{{ntf}Cp>qD(MuAB6L4-HmP+zO}1MkD@VGXE8k3=q6e?W9~(_Dp?N#;{yW99=5iu&UY*FR~R)riakkOu6BAlLd?2|@-k~%g72NR3zb4zCp zmM2q^k~26S?OsDoBUE?QlWU5Vr(aPUF@o$yH!(5jcOF{TAGZ^86$=tO6R_|NFmU9? zyQB{31;aF4xd+}+_M@eMrTqJ`IGu_o{C?NzXK0dM-$=SO(u<~L>3oai9u=SF>>3Fv zK#cj=E=Z!ebB3BL=0B&4&3p6E*s@Myz~FmQe;UAJ$MLAw`xYOT-i{Mp075SIc!JX} zh{67Wn=16W4cXZW`ng%XF3lN=S$pjCyBK1ypOZE7lHoYCNpaJDe%4fL>_!h`pV5?w zNPgS%t_)qh4x2;{eK2JuQ#`$x;?W5fSGLq+(qO(*-x`&9?Z8)pFA=oX!mTu197j z)}QFl+|f3`!rg2d#n3RCLI*-Vy}h$bn`Y=wrfX7tt}p5>`Q>^qF)#m6mJzSvGC^UI zxL5`Q6ih_OXRv7Isb(GJvTX)IN)k_qi^|j{yh~<(Rl57y+IkJ^d_Nzt9SL@S^iCP3 z5%{rYy@E(O78_F($!qQ7@Qr>a)f-(+*WpzP!_PqLu_1IMycts0sxNSfdzxeN*vEoa z+neERax>1}3@$oj(&rQToEKT>U5p6@rFOZl)ky7!s5GPfJ=~2k|Hw7lV_*;(Xjngkp!y>FEsO@C72O)%SQAlPw- zwqVexXZYvgEwr1Ia}y{@b9Sm96AUPCPN`?hez=0^ZI^ji53L7ix9kWbi3e^znt-Cc=M!4oBRsNG``7JI$YJGr2Y?7`$3 zSw%7@T2P1BO{>Merj!m1clz!{G~zX_+0zbaW4MHY!9-7gEUlgJ-W{IVY~&=B7#5QB zXUw$k2O}TWD?swv6nOJN061=#_Tiezptvkt&xP($7jo3)yDc^ucqMUk zOH(7((J^4hjx2u)5^w!|;$)?VjD-ijHzbqBOm|9Nh&6isl6?mMltFj)B=M@-OpZrn60W zL26GQK6RK(Z-)Xpll`sEqYR&B{O>a7Rv~qApE~D^s7nw{hu5A+=i6KJ)!cjIIM3e& z9V+d0aN)&`Qv6O!>xT8N=45tHWllI}ha)Q$|hHO;8Z|PTsSzA(#>1w1wlXWu95_TG&C? zy4I1xTtR%h*tRAHLRRXWRaWX64BxHQ8*h5s)g>7Di6xg>KUY z#?UfhLY=4!;^15_YnPVTZ-|=}>e5;w!`Vx}qBXOPjvA7P;B> zRe&8%oq7{DB6XiWX85KW@maw*n7aOuX<$}dZNRV9jg1-A7o8%A5cM{_80o-dR!?0N)v-8@&0k}2>dM%L(W}#rof^o} zBx@dgANswfYcLkaG($JfHn^U)EEFbn;yc~n={=y|fK;&CMv};InTN}BZ&JKmrr2U` zeiw=jCwgLfU>b@E))o8_i&XJ<=Tw_y6>ujnhjCYv(htdo3@>MHn$H^)*3~&}1>FzS z1Wo};0eeR5w$S)SEhgaea0APEF>JPxc-h?mzlb26e{zAElpZwA0poeTcjKlSxuX>2 zknyz;?uhA1;yi_y;OyB>=mlPv3B|U^p8$c2mPKJ#mTrwg%D^|Zj1C`M@9}F+I(Frn zXtX`*%>^KRZ@EnuCV~Opd)t48CIkNh(U)N=&Yy{MKGjIyST(NSl<+r_84SX)3@?OtH)^yAt}OC@|dgS~qkd(OttB zz}*og_xA72F3GI^@_bCaN;2=FdSFcui$~hu&FN(a8e2}&pojA2Bak=UMS>xGlIsoH z9=qap+uE>rE^ZN-gLarZs>k>8svp-jHMxn%lro|-4xA^};6?x5noCw@QJFK%wiPFI zJR@C(TA{leUFikYBk-I4V6MAY8YtFZ5ltK3d2A5XIrB?BS@T(_Y7i~v=D&vzMo>j* zcw#wTeZ=NAj$EVS=IJ2~SH%69ZESH!1vAx2_6W1=&u~$aTLMyOsZugF^#>A&sU3f9c}0(2&a*dO zy15^0bsI!wS1{0%6y`mCOuj3c#a$r^EbW;c9Ci+^-vh1+fHZAJD@o#-Cqk$*R}UG}TD zpHD}g`D$13H0%UF|6gyvjnjCqw_Ul65q`^`wgcGVj|)cs;(*V*%2%al{6`suH=-HE z5Q{zc@@LOkzronuB?Ub5VM@l1itr={PV3`8A#c8|e~Ss^UmnrpN-hcI;c$W}!nry5 z%cb{EaK1}=??lq!_%|VyOnDw~x*SHIhtqN9!ckWgT$z7W-qm58cQrU^JDX?IaNJ|+ zL4=4}FyKP$s5ixUuXJHg3|heZD1@YY`%Ddg?ZTh9@QW_|tP5{-;SDal)P?(9_#!r? zli#;o_-)(_$oFSl_y9u6@p%_^xp12c8%(yGkWu3EhzR+75+UQWJdoE@T$uMKWSaQj zyYOc&`{Gtm#?ZR7Kc%2KcaN+qbH@S`rg)`h!V*zdwq zT)4)Co*m#ki|}PVCV$p9h7lJckawID5xp#MVWkU6k1P~ScA;laj88ZMKg8p>_4+=$ z8R1yeC&~en<5y@6>e#YMgg9|NyBr~kaOy({Q4_Z@1kZ8zAw>P%$`J6emm-9nZes{O z<1TV|hJcUVj1WR>V~9Uq6ff}^0zP&YLL7YE#t`z3o8j;b0UwK_A7&grD0dD+D8@T& zd6j$j4ro5ccV_+z_>aHm@t&CxyZ71d{h)h43G_Vb zSgFxXAKyX$w8H~M`TalYH}l=)&&N;6!Ban(^q9BH>ZmofqYJyLk&Pk_qXlA zh6Q^a^NbK_wM6pA1vquGdv_G8Dye>4podEyqdOanH*A=bz_j8>qlx!Ae-7z4rq{IK lc2>-Qb1>%ZENN}Srcg}fH!_3oX-@tyvME^it?(It{|9usH30ws diff --git a/register/samplers/libsampler.cpp b/register/samplers/libsampler.cpp deleted file mode 100644 index fec46c6..0000000 --- a/register/samplers/libsampler.cpp +++ /dev/null @@ -1,285 +0,0 @@ -#include -#include "ndarray.h" -#include "math.h" - -using namespace std; - -extern "C" { - -#include - -/* - -A mapping function which adjusts coordinates outside the domain in the following way: - - 'n' nearest : sampler the nearest valid coordinate. - - 'r' reflect : reflect the coordinate into the boundary. - - 'w' wrap : wrap the coordinate around the appropriate axis. -*/ - -int map(int *coords, int rows, int cols, char mode) - { - switch (mode) - { - - case 'c': /* constant */ - - if (coords[0] < 0) return 0; - if (coords[1] < 0) return 0; - - if (coords[0] == rows) return 0; - if (coords[1] == cols) return 0; - - if (coords[0] > rows) return 0; - if (coords[1] > cols) return 0; - - break; - - case 'n': /* nearest */ - - if (coords[0] < 0) coords[0] = 0; - if (coords[1] < 0) coords[1] = 0; - - if (coords[0] == rows) coords[0] = rows-1; - if (coords[1] == cols) coords[1] = cols-1; - - if (coords[0] > rows) coords[0] = rows; - if (coords[1] > cols) coords[1] = cols; - - break; - - case 'r': /* reflect */ - - if (coords[0] < 0) coords[0] = fmod(-coords[0],rows); - if (coords[1] < 0) coords[1] = fmod(-coords[1],cols); - - if (coords[0] == rows) coords[0] = rows-1; - if (coords[1] == cols) coords[1] = cols-1; - - if (coords[0] > rows) coords[0] = rows; - if (coords[1] > cols) coords[1] = cols; - - break; - - case 'w': /* wrap */ - - if (coords[0] < 0) coords[0] = rows - fmod(-coords[0],rows); - if (coords[1] < 0) coords[1] = cols - fmod(-coords[1],cols); - - if (coords[0] == rows) coords[0] = 0; - if (coords[1] == cols) coords[1] = 0; - - if (coords[0] > rows) coords[0] = fmod(coords[0],rows); - if (coords[1] > cols) coords[1] = fmod(coords[1],cols); - - break; - } - - return 1; - } - -/* - nearest neighbour interpolator. - */ - -int nearest(numpyArray array0, - numpyArray array1, - numpyArray array2, - char mode, - double cvalue - ) -{ - Ndarray warp(array0); - Ndarray image(array1); - Ndarray result(array2); - - int coords[2] = {0, 0}; - - int rows = image.getShape(0); - int cols = image.getShape(1); - - for (int i = 0; i < warp.getShape(1); i++) - { - for (int j = 0; j < warp.getShape(2); j++) - { - coords[0] = (int)warp[0][i][j]; - coords[1] = (int)warp[1][i][j]; - - if ( not map(coords, rows, cols, mode) ) - { - result[i][j] = cvalue; - } - else - { - result[i][j] = image[coords[0]][coords[1]]; - } - } - } - - return 0; -} - - -/* - bilinear interpolator - */ - -int bilinear(numpyArray array0, - numpyArray array1, - numpyArray array2, - char mode, - double cvalue - ) -{ - Ndarray warp(array0); - Ndarray image(array1); - Ndarray result(array2); - - double di = 0.0; - double dj = 0.0; - - double fi = 0; - double fj = 0; - - double w0 = 0.0; - double w1 = 0.0; - double w2 = 0.0; - double w3 = 0.0; - - int tl[2] = {0, 0}; - int tr[2] = {0, 0}; - int ll[2] = {0, 0}; - int lr[2] = {0, 0}; - - int rows = image.getShape(0); - int cols = image.getShape(1); - - for (int i = 0; i < warp.getShape(1); i++) - { - for (int j = 0; j < warp.getShape(2); j++) - { - /* Floating point coordinates */ - fi = warp[0][i][j]; - fj = warp[1][i][j]; - - /* Integer component */ - di = (double)((int)(warp[0][i][j])); - dj = (double)((int)(warp[1][i][j])); - - /* Defined sampling coordinates */ - - tl[0] = (int)fi; - tl[1] = (int)fj; - - tr[0] = tl[0]; - tr[1] = tl[1] + 1; - - ll[0] = tl[0] + 1; - ll[1] = tl[1]; - - lr[0] = tl[0] + 1; - lr[1] = tl[1] + 1; - - w0 = 0.0; - if ( map(tl, rows, cols, mode) ) - w0 = ((dj+1-fj)*(di+1-fi))*image[tl[0]][tl[1]]; - - w1 = 0.0; - if ( map(tr, rows, cols, mode) ) - w1 = ((fj-dj)*(di+1-fi))*image[tr[0]][tr[1]]; - - w2 = 0.0; - if ( map(ll, rows, cols, mode) ) - w2 = ((dj+1-fj)*(fi-di))*image[ll[0]][ll[1]]; - - w3 = 0.0; - if ( map(lr, rows, cols, mode) ) - w3 = ((fj-dj)*(fi-di))*image[lr[0]][lr[1]]; - - result[i][j] = w0 + w1 + w2 + w3; - - } - } - - return 0; -} - - -int cubicConvolution(numpyArray array0, - numpyArray array1, - numpyArray array2 - ) -{ - Ndarray warp(array0); - Ndarray image(array1); - Ndarray result(array2); - - int di = 0; - int dj = 0; - - int rows = image.getShape(0); - int cols = image.getShape(1); - - double xShift; - double yShift; - double xArray0; - double xArray1; - double xArray2; - double xArray3; - double yArray0; - double yArray1; - double yArray2; - double yArray3; - double c0; - double c1; - double c2; - double c3; - - for (int i = 0; i < image.getShape(0); i++) - { - for (int j = 0; j < image.getShape(1); j++) - { - di = (int)floor(warp[0][i][j]); - dj = (int)floor(warp[1][i][j]); - - if ( ( di < rows-2 && di >= 2 ) && - ( dj < cols-2 && dj >= 2 ) ) - { - xShift = warp[1][i][j] - dj; - yShift = warp[0][i][j] - di; - xArray0 = -(1/2.0)*pow(xShift, 3) + pow(xShift, 2) - (1/2.0)*xShift; - xArray1 = (3/2.0)*pow(xShift, 3) - (5/2.0)*pow(xShift, 2) + 1; - xArray2 = -(3/2.0)*pow(xShift, 3) + 2*pow(xShift, 2) + (1/2.0)*xShift; - xArray3 = (1/2.0)*pow(xShift, 3) - (1/2.0)*pow(xShift, 2); - yArray0 = -(1/2.0)*pow(yShift, 3) + pow(yShift, 2) - (1/2.0)*yShift; - yArray1 = (3/2.0)*pow(yShift, 3) - (5/2.0)*pow(yShift, 2) + 1; - yArray2 = -(3/2.0)*pow(yShift, 3) + 2*pow(yShift, 2) + (1/2.0)*yShift; - yArray3 = (1/2.0)*pow(yShift, 3) - (1/2.0)*pow(yShift, 2); - c0 = xArray0 * image[di-1][dj-1] + xArray1 * image[di-1][dj+0] + xArray2 * image[di-1][dj+1] + xArray3 * image[di-1][dj+2]; - c1 = xArray0 * image[di+0][dj-1] + xArray1 * image[di+0][dj+0] + xArray2 * image[di+0][dj+1] + xArray3 * image[di+0][dj+2]; - c2 = xArray0 * image[di+1][dj-1] + xArray1 * image[di+1][dj+0] + xArray2 * image[di+1][dj+1] + xArray3 * image[di+1][dj+2]; - c3 = xArray0 * image[di+2][dj-1] + xArray1 * image[di+2][dj+0] + xArray2 * image[di+2][dj+1] + xArray3 * image[di+2][dj+2]; - result[i][j] = c0 * yArray0 + c1 * yArray1 + c2 * yArray2 + c3 * yArray3; - } - else - { - result[i][j] = 0.0; - } - } - } - - return 0; -} - -static PyMethodDef sampler_methods[] = { - {NULL, NULL} -}; - -void initlibsampler() -{ - (void) Py_InitModule("libsampler", sampler_methods); -} - -} // end extern "C" diff --git a/register/samplers/ndarray.h b/register/samplers/ndarray.h deleted file mode 100644 index 0ecf36a..0000000 --- a/register/samplers/ndarray.h +++ /dev/null @@ -1,206 +0,0 @@ -//========================================================== -// ndarray.h: C++ interface for easy access to numpy arrays -// -// J. De Ridder -//========================================================== - - -#ifndef NDARRAY_H -#define NDARRAY_H - - - -// The C-struct to retrieve the ctypes structure. -// Note: Order of struct members must be the same as in Python. -// Member names are not recognized! - - -template -struct numpyArray -{ - T *data; - long *shape; - long *strides; -}; - - - - -// Traits need to used because the return type of the []-operator can be -// a subarray or an element, depending whether all the axes are exhausted. - -template class Ndarray; // forward declaration - - -template -struct getItemTraits -{ - typedef Ndarray returnType; -}; - - -template -struct getItemTraits -{ - typedef datatype& returnType; -}; - - - -// Ndarray definition - -template -class Ndarray -{ - private: - datatype *data; - long *shape; - long *strides; - - public: - Ndarray(datatype *data, long *shape, long *strides); - Ndarray(const Ndarray& array); - Ndarray(const numpyArray& array); - long getShape(const int axis); - typename getItemTraits::returnType operator[](unsigned long i); -}; - - -// Ndarray constructor - -template -Ndarray::Ndarray(datatype *data, long *shape, long *strides) -{ - this->data = data; - this->shape = shape; - this->strides = strides; -} - - -// Ndarray copy constructor - -template -Ndarray::Ndarray(const Ndarray& array) -{ - this->data = array.data; - this->shape = array.shape; - this->strides = array.strides; -} - - -// Ndarray constructor from ctypes structure - -template -Ndarray::Ndarray(const numpyArray& array) -{ - this->data = array.data; - this->shape = array.shape; - this->strides = array.strides; -} - - -// Ndarray method to get length of given axis - -template -long Ndarray::getShape(const int axis) -{ - return this->shape[axis]; -} - - - -// Ndarray overloaded []-operator. -// The [i][j][k] selection is recursively replaced by i*strides[0]+j*strides[1]+k*strides[2] -// at compile time, using template meta-programming. If the axes are not exhausted, return -// a subarray, else return an element. - -template -typename getItemTraits::returnType -Ndarray::operator[](unsigned long i) -{ - return Ndarray(&this->data[i*this->strides[0]], &this->shape[1], &this->strides[1]); -} - - - -// Template partial specialisation of Ndarray. -// For 1D Ndarrays, the [] operator should return an element, not a subarray, so it needs -// to be special-cased. In principle only the operator[] method should be specialised, but -// for some reason my gcc version seems to require that then the entire class with all its -// methods are specialised. - -template -class Ndarray -{ - private: - datatype *data; - long *shape; - long *strides; - - public: - Ndarray(datatype *data, long *shape, long *strides); - Ndarray(const Ndarray& array); - Ndarray(const numpyArray& array); - long getShape(const int axis); - typename getItemTraits::returnType operator[](unsigned long i); -}; - - -// Ndarray partial specialised constructor - -template -Ndarray::Ndarray(datatype *data, long *shape, long *strides) -{ - this->data = data; - this->shape = shape; - this->strides = strides; -} - - - -// Ndarray partially specialised copy constructor - -template -Ndarray::Ndarray(const Ndarray& array) -{ - this->data = array.data; - this->shape = array.shape; - this->strides = array.strides; -} - - - -// Ndarray partially specialised constructor from ctypes structure - -template -Ndarray::Ndarray(const numpyArray& array) -{ - this->data = array.data; - this->shape = array.shape; - this->strides = array.strides; -} - - - -// Ndarray method to get length of given axis - -template -long Ndarray::getShape(const int axis) -{ - return this->shape[axis]; -} - - - -// Partial specialised [] operator: for 1D arrays, return an element rather than a subarray - -template -typename getItemTraits::returnType -Ndarray::operator[](unsigned long i) -{ - return this->data[i*this->strides[0]]; -} - - - -#endif diff --git a/register/samplers/numpyctypes.py b/register/samplers/numpyctypes.py deleted file mode 100644 index 5209a8f..0000000 --- a/register/samplers/numpyctypes.py +++ /dev/null @@ -1,109 +0,0 @@ -import numpy as N -import ctypes as C - - -ctypesDict = {'d' : C.c_double, - 'b' : C.c_char, - 'h' : C.c_short, - 'i' : C.c_int, - 'l' : C.c_long, - 'q' : C.c_longlong, - 'B' : C.c_ubyte, - 'H' : C.c_ushort, - 'I' : C.c_uint, - 'L' : C.c_ulong, - 'Q' : C.c_ulonglong} - - -def c_ndarray(a, dtype = None, ndim = None, shape = None, requirements = None): - - """ - PURPOSE: Given an array, return a ctypes structure containing the - arrays info (data, shape, strides, ndim). A check is made to ensure - that the array has the specified dtype and requirements. - - INPUT: a: an array: something that is or can be converted to a numpy array - dtype: the required dtype of the array, convert if it doesn't match - ndim: integer: the required number of axes of the array - shape: tuple of integers: required shape of the array - requirements: list of requirements: (E)nsurearray, (F)ortran, (F)_contiguous, - (C)ontiguous, (C)_contiguous. Convert if it doesn't match. - - OUTPUT: ctypes structure with the fields: - . data: pointer to the data : the type is determined with the dtype of the - array, and with ctypesDict. - . shape: pointer to long array : size of each of the dimensions - . strides: pointer to long array : strides in elements (not bytes) - """ - - if not requirements: - - # Also allow derived classes of ndarray - array = N.asanyarray(a, dtype=dtype) - else: - - # Convert requirements to captial letter codes: - # (ensurearray' -> 'E'; 'aligned' -> 'A' - # 'fortran', 'f_contiguous', 'f' -> 'F' - # 'contiguous', 'c_contiguous', 'c' -> 'C') - - requirements = [x[0].upper() for x in requirements] - subok = (0 if 'E' in requirements else 1) - - # Make from 'a' an ndarray with the specified dtype, but don't copy the - # data (yet). This also ensures that the .flags attribute is present. - - array = N.array(a, dtype=dtype, copy=False, subok=subok) - - # See if copying all data is really necessary. - # Note: 'A' = (A)ny = only (F) it is was already (F) - - copychar = 'A' - if 'F' in requirements: - copychar = 'F' - elif 'C' in requirements: - copychar = 'C' - - for req in requirements: - if not array.flags[req]: - array = array.copy(copychar) - break - - - # If required, check the number of axes and the shape of the array - - if ndim is not None: - if array.ndim != ndim: - raise TypeError, "Array has wrong number of axes" - - if shape is not None: - if array.shape != shape: - raise TypeError, "Array has wrong shape" - - # Define a class that serves as interface of an ndarray to ctypes. - # Part of the type depends on the array's dtype. - - class ndarrayInterfaceToCtypes(C.Structure): - pass - - typechar = array.dtype.char - - if typechar in ctypesDict: - ndarrayInterfaceToCtypes._fields_ = \ - [("data", C.POINTER(ctypesDict[typechar])), - ("shape" , C.POINTER(C.c_long)), - ("strides", C.POINTER(C.c_long))] - else: - raise TypeError, "dtype of input ndarray not supported" - - # Instantiate the interface class and attach the ndarray's internal info. - # Ctypes does automatic conversion between (c_long * #) arrays and POINTER(c_long). - - ndarrayInterface = ndarrayInterfaceToCtypes() - ndarrayInterface.data = array.ctypes.data_as(C.POINTER(ctypesDict[typechar])) - ndarrayInterface.shape = (C.c_long * array.ndim)(*array.shape) - ndarrayInterface.strides = (C.c_long * array.ndim)(*array.strides) - for n in range(array.ndim): - ndarrayInterface.strides[n] /= array.dtype.itemsize - - return ndarrayInterface diff --git a/register/samplers/sampler.py b/register/samplers/sampler.py deleted file mode 100644 index 4068269..0000000 --- a/register/samplers/sampler.py +++ /dev/null @@ -1,266 +0,0 @@ -""" A collection of samplers""" - -import numpy as np -import scipy.ndimage as nd - -from numpy.ctypeslib import load_library -from numpyctypes import c_ndarray -import ctypes - -libsampler = load_library('libsampler', __file__) - -# Configuration for the extrapolation mode and fill value. -EXTRAPOLATION_MODE = 'c' -EXTRAPOLATION_CVALUE = 0.0 - - -class Sampler(object): - """ - Abstract sampler - - Attributes - ---------- - METRIC : string - The type of similarity sampler being used. - DESCRIPTION : string - A meaningful description of the sampler used, with references where - appropriate. - """ - - METHOD=None - DESCRIPTION=None - - def __init__(self, coordinates): - - self.coordinates = coordinates - - def f(self, array, warp): - """ - A sampling function, responsible for returning a sampled set of values - from the given array. - - Parameters - ---------- - array: nd-array - Input array for sampling. - warp: nd-array - Deformation coordinates. - - Returns - ------- - sample: nd-array - Sampled array data. - """ - - if self.coordinates is None: - raise ValueError('Appropriately defined coordinates not provided.') - - i = self.coordinates.tensor[0] + warp[0] - j = self.coordinates.tensor[1] + warp[1] - - packedCoords = (i.reshape(1,i.size), - j.reshape(1,j.size)) - - return self.sample(array, np.vstack(packedCoords)) - - def sample(self, array, coords): - """ - The sampling function - provided by the specialized samplers. - """ - return None - - def __str__(self): - return 'Method: {0} \n {1}'.format( - self.METHOD, - self.DESCRIPTION - ) - - -class Nearest(Sampler): - - METHOD='Nearest Neighbour (NN)' - - DESCRIPTION=""" - Given coordinate in the array nearest neighbour sampling simply rounds - coordinates points: - f(I; i,j) = I( round(i), round(j)) - """ - - def __init__(self, coordinates): - Sampler.__init__(self, coordinates) - - - def f(self, array, warp): - """ - A sampling function, responsible for returning a sampled set of values - from the given array. - - Parameters - ---------- - array: nd-array - Input array for sampling. - warp: nd-array - Deformation coordinates. - - Returns - ------- - sample: nd-array - Sampled array data. - """ - - if self.coordinates is None: - raise ValueError('Appropriately defined coordinates not provided.') - - result = np.zeros_like(warp[0]) - - arg0 = c_ndarray(warp, dtype=np.double, ndim=3) - arg1 = c_ndarray(array, dtype=np.double, ndim=2) - arg2 = c_ndarray(result, dtype=np.double, ndim=2) - - libsampler.nearest( - arg0, - arg1, - arg2, - ctypes.c_char(EXTRAPOLATION_MODE[0]), - ctypes.c_double(EXTRAPOLATION_CVALUE) - ) - - return result.flatten() - - -class Bilinear(Sampler): - - METHOD='Bilinear (BL)' - - DESCRIPTION=""" - Given a coordinate in the array a linear interpolation is performed - between 4 (2x2) nearest values. - """ - - def __init__(self, coordinates): - Sampler.__init__(self, coordinates) - - def f(self, array, warp): - """ - A sampling function, responsible for returning a sampled set of values - from the given array. - - Parameters - ---------- - array: nd-array - Input array for sampling. - warp: nd-array - Deformation coordinates. - - Returns - ------- - sample: nd-array - Sampled array data. - """ - - if self.coordinates is None: - raise ValueError('Appropriately defined coordinates not provided.') - - result = np.zeros_like(warp[0]) - - arg0 = c_ndarray(warp, dtype=np.double, ndim=3) - arg1 = c_ndarray(array, dtype=np.double, ndim=2) - arg2 = c_ndarray(result, dtype=np.double, ndim=2) - - libsampler.bilinear( - arg0, - arg1, - arg2, - ctypes.c_char(EXTRAPOLATION_MODE[0]), - ctypes.c_double(EXTRAPOLATION_CVALUE) - ) - - return result.flatten() - - -class CubicConvolution(Sampler): - - METHOD='Cubic Convolution (CC)' - - DESCRIPTION=""" - Given a coordinate in the array cubic convolution interpolates between - 16 (4x4) nearest values. - """ - - def __init__(self, coordinates): - Sampler.__init__(self, coordinates) - - - def f(self, array, warp): - """ - A sampling function, responsible for returning a sampled set of values - from the given array. - - Parameters - ---------- - array: nd-array - Input array for sampling. - warp: nd-array - Deformation coordinates. - - Returns - ------- - sample: nd-array - Sampled array data. - """ - - if self.coordinates is None: - raise ValueError('Appropriately defined coordinates not provided.') - - result = np.zeros_like(array) - - arg0 = c_ndarray(warp, dtype=np.double, ndim=3) - arg1 = c_ndarray(array, dtype=np.double, ndim=2) - arg2 = c_ndarray(result, dtype=np.double, ndim=2) - - libsampler.cubicConvolution(arg0, arg1, arg2) - - return result.flatten() - - -class Spline(Sampler): - - METHOD='nd-image spline sampler (SR)' - - DESCRIPTION=""" - Refer to the documentation for the ndimage map_coordinates function. - - http://docs.scipy.org/doc/scipy/reference/generated/ - scipy.ndimage.interpolation.map_coordinates.html - """ - - def __init__(self, coordinates): - Sampler.__init__(self, coordinates) - - def f(self, array, warp): - """ - A sampling function, responsible for returning a sampled set of values - from the given array. - - Parameters - ---------- - array: nd-array - Input array for sampling. - warp: nd-array - Deformation coordinates. - - Returns - ------- - sample: nd-array - Sampled array data. - """ - - if self.coordinates is None: - raise ValueError('Appropriately defined coordinates not provided.') - - return nd.map_coordinates( - array, - warp, - order=2, - mode='nearest' - ).flatten() diff --git a/register/samplers/setup.py b/register/samplers/setup.py deleted file mode 100644 index 36f6f85..0000000 --- a/register/samplers/setup.py +++ /dev/null @@ -1,24 +0,0 @@ -#!/usr/bin/env python - -import os - -base_path = os.path.abspath(os.path.dirname(__file__)) - -def configuration(parent_package='', top_path=None): - from numpy.distutils.misc_util import \ - Configuration, get_numpy_include_dirs - - config = Configuration('samplers', parent_package, top_path) - - config.add_extension('libsampler', sources=['libsampler.cpp'], - include_dirs=[get_numpy_include_dirs()]) - - return config - - -if __name__ == '__main__': - from numpy.distutils.core import setup - - setup( - **(configuration(top_path='').todict()) - ) diff --git a/register/setup.py b/register/setup.py deleted file mode 100644 index c558394..0000000 --- a/register/setup.py +++ /dev/null @@ -1,21 +0,0 @@ -import os - -def configuration(parent_package='', top_path=None): - - from numpy.distutils.misc_util import Configuration - - config = Configuration('register', parent_package, top_path) - - config.add_subpackage('models') - config.add_subpackage('metrics') - config.add_subpackage('samplers') - config.add_subpackage('visualize') - config.add_subpackage('features') - - return config - -if __name__ == "__main__": - from numpy.distutils.core import setup - - config = configuration(top_path='').todict() - setup(**config) diff --git a/register/visualize/__init__.py b/register/visualize/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/visualize/dialog.py b/register/visualize/dialog.py deleted file mode 100644 index 360710d..0000000 --- a/register/visualize/dialog.py +++ /dev/null @@ -1,158 +0,0 @@ -""" A custom dialog for introspection of search paths.""" - -import numpy as np -from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas -from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar -from matplotlib.figure import Figure - -from PyQt4 import QtCore, QtGui - -#=============================================================================== -# Custom dialog generated by qt-designer. -#=============================================================================== - -try: - _fromUtf8 = QtCore.QString.fromUtf8 -except AttributeError: - _fromUtf8 = lambda s: s - -class Ui_Dialog(object): - def setupUi(self, Dialog): - Dialog.setObjectName(_fromUtf8("Dialog")) - Dialog.resize(1233, 624) - Dialog.setWindowTitle(QtGui.QApplication.translate("Dialog", "Dialog", None, QtGui.QApplication.UnicodeUTF8)) - self.horizontalLayout = QtGui.QHBoxLayout(Dialog) - self.horizontalLayout.setObjectName(_fromUtf8("horizontalLayout")) - self.splitter = QtGui.QSplitter(Dialog) - sizePolicy = QtGui.QSizePolicy(QtGui.QSizePolicy.Expanding, QtGui.QSizePolicy.Preferred) - sizePolicy.setHorizontalStretch(0) - sizePolicy.setVerticalStretch(0) - sizePolicy.setHeightForWidth(self.splitter.sizePolicy().hasHeightForWidth()) - self.splitter.setSizePolicy(sizePolicy) - self.splitter.setMinimumSize(QtCore.QSize(400, 0)) - self.splitter.setOrientation(QtCore.Qt.Horizontal) - self.splitter.setObjectName(_fromUtf8("splitter")) - self.widget = MyDynamicMplCanvas(self.splitter) - sizePolicy = QtGui.QSizePolicy(QtGui.QSizePolicy.Maximum, QtGui.QSizePolicy.Preferred) - sizePolicy.setHorizontalStretch(0) - sizePolicy.setVerticalStretch(0) - sizePolicy.setHeightForWidth(self.widget.sizePolicy().hasHeightForWidth()) - self.widget.setSizePolicy(sizePolicy) - self.widget.setMinimumSize(QtCore.QSize(1000, 600)) - self.widget.setObjectName(_fromUtf8("widget")) - - self.layoutWidget = QtGui.QWidget(self.splitter) - self.layoutWidget.setObjectName(_fromUtf8("layoutWidget")) - self.verticalLayout_2 = QtGui.QVBoxLayout(self.layoutWidget) - self.verticalLayout_2.setSpacing(10) - self.verticalLayout_2.setContentsMargins(-1, 6, -1, -1) - self.verticalLayout_2.setObjectName(_fromUtf8("verticalLayout_2")) - self.label = QtGui.QLabel(self.layoutWidget) - self.label.setMaximumSize(QtCore.QSize(200, 16777215)) - self.label.setText(QtGui.QApplication.translate("Dialog", "Fittings", None, QtGui.QApplication.UnicodeUTF8)) - self.label.setObjectName(_fromUtf8("label")) - self.verticalLayout_2.addWidget(self.label) - self.listWidget = QtGui.QListWidget(self.layoutWidget) - self.listWidget.setObjectName(_fromUtf8("listWidget")) - self.verticalLayout_2.addWidget(self.listWidget) - self.horizontalLayout.addWidget(self.splitter) - self.verticalLayout_2.addWidget(NavigationToolbar(self.widget, Dialog)) - - self.listWidget.connect( - self.listWidget, - QtCore.SIGNAL("currentItemChanged (QListWidgetItem *,QListWidgetItem *)"), - self.itemChangedSlot) - - self.retranslateUi(Dialog) - QtCore.QMetaObject.connectSlotsByName(Dialog) - - def itemChangedSlot(self, item, item2): - self.widget.formPlot(self.searchMap[item]) - - def updateList(self, search): - """ - Takes the search datastructure and makes a list of it. - """ - - self.searchMap = {} - for index, step in enumerate(search): - item = QtGui.QListWidgetItem(self.listWidget) - item.setText( - "#{0}, e: {1}".format( - index, - step.error - ) - ) - - if not step.decreasing: - item.setBackgroundColor(QtGui.QColor("red")) - self.searchMap[item] = step - - def retranslateUi(self, Dialog): - pass - -############################################################################### -# MPL widget -############################################################################### - -class MyMplCanvas(FigureCanvas): - """Ultimately, this is a QWidget (as well as a FigureCanvasAgg, etc.).""" - def __init__(self, parent=None, step=None): - self.fig = Figure() - - self.ax1 = self.fig.add_subplot(1,3,1) - self.ax2 = self.fig.add_subplot(1,3,2) - self.ax3 = self.fig.add_subplot(1,3,3)#, sharex=self.ax1, sharey=self.ax1) - - if step is not None: - self.formPlot(step) - - FigureCanvas.__init__(self, self.fig) - self.setParent(parent) - - FigureCanvas.setSizePolicy( - self, - QtGui.QSizePolicy.Expanding, - QtGui.QSizePolicy.Expanding) - - FigureCanvas.updateGeometry(self) - -class MyDynamicMplCanvas(MyMplCanvas): - """A canvas that updates itself every second with a new plot.""" - def __init__(self, *args, **kwargs): - MyMplCanvas.__init__(self, *args, **kwargs) - - def formPlot(self, step): - """ - A graphical representation of an optimization step. - - Parameters - ---------- - step: Register.optstep - A step in the registration. - """ - self.ax1.cla() - self.img1 = self.ax1.imshow( - step.image, - interpolation='nearest', - cmap='gray' - ) - self.ax1.axis('image') - - self.ax2.cla() - self.img2 = self.ax2.imshow( - step.warpedImage, - interpolation='nearest', - cmap='gray' - ) - self.ax2.axis('image') - - self.ax3.cla() - self.img3 = self.ax3.imshow( - step.template, - interpolation='nearest', - cmap='gray' - ) - self.ax2.axis('image') - - self.draw() diff --git a/register/visualize/dialog.ui b/register/visualize/dialog.ui deleted file mode 100644 index 2920c07..0000000 --- a/register/visualize/dialog.ui +++ /dev/null @@ -1,80 +0,0 @@ - - - Dialog - - - - 0 - 0 - 1233 - 624 - - - - Dialog - - - - - - - 0 - 0 - - - - - 400 - 0 - - - - Qt::Horizontal - - - - - 0 - 0 - - - - - 1000 - 600 - - - - - - - 10 - - - 6 - - - - - - 200 - 16777215 - - - - Fittings - - - - - - - - - - - - - - - diff --git a/register/visualize/plot.py b/register/visualize/plot.py deleted file mode 100644 index 7a916dc..0000000 --- a/register/visualize/plot.py +++ /dev/null @@ -1,225 +0,0 @@ -''' (Debug utility) Defines a set of plotting callback functions ''' - -import sys -import matplotlib.pyplot as plt - -IMAGE_ORIGIN=None -IMAGE_COLORMAP='gray' -IMAGE_VMIN=None -IMAGE_VMAX=None - -#=============================================================================== -# Plot configuration -#=============================================================================== - -params = { - 'axes.labelsize': 10, - 'axes.titlesize': 10, - 'figure.titlesize': 12, - 'font.size': 10, - 'font.weight':'normal', - 'text.fontsize': 10, - 'axes.fontsize': 10, - 'legend.fontsize': 11, - 'xtick.labelsize': 8, - 'ytick.labelsize': 8, - 'figure.figsize': (12,6), - 'figure.facecolor': 'w', - } - -plt.rcParams.update(params); - -#=============================================================================== -# Debug tools that make use of PyQt4. -#=============================================================================== - -def searchInspector(search): - """ - A Qt based GUI to inspect the output of search algorithms. - - Parameters - ---------- - search: list of Register.optstepnd-array - The output of a registration. - """ - - try: - from PyQt4.QtGui import QApplication, QDialog - from dialog import Ui_Dialog - except Exception: - print "Missing a required library - please install pyQt4." - return - - app = QApplication(sys.argv) - window = QDialog() - ui = Ui_Dialog() - ui.setupUi(window) - ui.updateList(search) - window.show() - app.exec_() - -#=============================================================================== -# Code below needs to be cleaned up. -#=============================================================================== - -plt.ion() - -def show(): - - plt.ioff() - plt.show() - -def coordPlt(grid, buffer=10, step=5): - """ - Plot the grid coordinates. - """ - plt.cla() - - plt.plot(grid[1][0::step, 0::step], - grid[0][0::step, 0::step], - '.-b' ) - - plt.plot(grid[1][0::step, 0::step].T, - grid[0][0::step, 0::step].T, - '.-b' ) - - plt.axis( [ grid[1].max() + buffer, - grid[1].min() - buffer, - grid[0].max() + buffer, - grid[0].min() - buffer], - ) - plt.axis('off') - plt.grid() - -def featurePlt(features): - - for id, point in features['points'].items(): - plt.plot(point[0], point[1], 'or') - - -def boundPlt(grid): - - xmin = grid[1].min() - ymin = grid[0].min() - - xmax = grid[1].max() - ymax = grid[0].max() - - plt.hlines([ymin,ymax], xmin, xmax, colors='g') - plt.vlines([xmin, xmax], ymin, ymax, colors='g') - -def warpPlot(grid, warp, _warp): - - plt.subplot(1,2,1) - coordPlt(warp) - boundPlt(grid) - - plt.subplot(1,2,2) - coordPlt(_warp) - boundPlt(grid) - - -def featurePlot(image, template=None, warpedImage=None): - - plt.subplot(1,4,1) - plt.title('I') - plt.imshow(image.data, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP, - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - featurePlt(image.features) - - if not template is None: - plt.subplot(1,3,2) - plt.title('T') - plt.imshow(template.data, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP, - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - featurePlt(template.features) - - if not warpedImage is None: - plt.subplot(1,3,3) - plt.title('W(I;p)') - plt.imshow(warpedImage, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP, - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - featurePlt(template.features) - -def featurePlotSingle(image): - plt.title('I') - plt.imshow(image.data, - cmap=IMAGE_COLORMAP, - origin='lower', - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - featurePlt(image.features) - - -def gridPlot(image, template, warpedImage, grid, warp, title): - - plt.subplot(2,3,1) - plt.title('I') - plt.imshow(image, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP, - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - - plt.subplot(2,3,2) - plt.title('T') - plt.imshow(template, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP, - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - - plt.subplot(2,3,3) - plt.title('W(I;p)') - plt.imshow(warpedImage, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP, - vmin=IMAGE_VMIN, - vmax=IMAGE_VMAX - ) - plt.axis('off') - - plt.subplot(2,3,4) - plt.title('I-T') - plt.imshow(template - image, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP - ) - plt.axis('off') - - plt.subplot(2,3,5) - plt.axis('off') - coordPlt(warp) - boundPlt(grid) - plt.title('W(x;p)') - - plt.subplot(2,3,6) - plt.title('W(I;p) - T {0}'.format(title)) - plt.imshow(template - warpedImage, - origin=IMAGE_ORIGIN, - cmap=IMAGE_COLORMAP - ) - plt.axis('off') - - plt.draw() diff --git a/setup.py b/setup.py index 5c60f46..560b141 100644 --- a/setup.py +++ b/setup.py @@ -4,8 +4,8 @@ """ -DISTNAME = 'python-register' -DESCRIPTION = 'Image registration toolbox for SciPy' +DISTNAME = 'imreg' +DESCRIPTION = '(imreg) Image registration' LONG_DESCRIPTION = descr MAINTAINER = 'Nathan Faggian' MAINTAINER_EMAIL = 'nathan.faggian@gmail.com' @@ -34,8 +34,7 @@ def configuration(parent_package='', top_path=None): delegate_options_to_subpackages=True, quiet=True) - config.add_subpackage('register') -# config.add_subpackage(DISTNAME) + config.add_subpackage('imreg') return config diff --git a/tests/test_detectors.py b/tests/test_detectors.py index b6aef18..2c8a042 100644 --- a/tests/test_detectors.py +++ b/tests/test_detectors.py @@ -1,6 +1,6 @@ import os import matplotlib.pyplot as plt -import register.features.detector as detector +import imreg.features.detector as detector def test_haardetector(): """ diff --git a/tests/test_haar2d.py b/tests/test_haar2d.py index f2e5fa2..60eda9a 100644 --- a/tests/test_haar2d.py +++ b/tests/test_haar2d.py @@ -1,5 +1,5 @@ import numpy as np -import register.features.haar2d as haar2d +import imreg.features.haar2d as haar2d def test_haar2d(): """ diff --git a/tests/test_models.py b/tests/test_models.py index a3eeff5..22d68f7 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -1,8 +1,8 @@ import numpy as np -import register.models.model as model -import register.samplers.sampler as sampler -import register.register as register +import imreg.models.model as model +import imreg.samplers.sampler as sampler +import imreg.register as register import scipy.misc as misc import scipy.ndimage as nd diff --git a/tests/test_register.py b/tests/test_register.py index 88ec02d..4dc6c21 100644 --- a/tests/test_register.py +++ b/tests/test_register.py @@ -3,11 +3,11 @@ import scipy.ndimage as nd import scipy.misc as misc -import register.models.model as model -import register.metrics.metric as metric -import register.samplers.sampler as sampler +import imreg.models.model as model +import imreg.metrics.metric as metric +import imreg.samplers.sampler as sampler -from register import register +from imreg import register def warp(image, p, model, sampler): """ diff --git a/tests/test_register_data.py b/tests/test_register_data.py index 1aca27e..fdc723c 100644 --- a/tests/test_register_data.py +++ b/tests/test_register_data.py @@ -1,6 +1,6 @@ import scipy.misc as misc -from register import register +from imreg import register def test_downsample(): """ diff --git a/tests/test_register_features.py b/tests/test_register_features.py index d322091..063b6f7 100644 --- a/tests/test_register_features.py +++ b/tests/test_register_features.py @@ -1,8 +1,8 @@ import numpy as np -from register.models import model -from register.samplers import sampler -from register import register +from imreg.models import model +from imreg.samplers import sampler +from imreg import register def test_register(): """ diff --git a/tests/test_samplers.py b/tests/test_samplers.py index 931eff6..b20c841 100644 --- a/tests/test_samplers.py +++ b/tests/test_samplers.py @@ -3,10 +3,10 @@ import numpy as np import scipy as sp -import register.models.model as model -import register.samplers.sampler as sampler +import imreg.models.model as model +import imreg.samplers.sampler as sampler -from register import register +from imreg import register def test_sampler():