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 diff --git a/LICENSE.txt b/LICENSE.txt index 08218b3..f4f65c3 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,4 +1,8 @@ -Copyright 2011, Nathan Faggian (nathan.faggian@gmail.com) +Copyright 2011, +Nathan Faggian, Australian Bureau of Meteorology. (nfaggian@bom.gov.au) +Stefan Van Der Walt. +Riaan Van Den Dool. (riaanvddool@gmail.com) + Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. diff --git a/README.rst b/README.rst index 633a9b0..5b25176 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 ~~~~~~~~~~~~~ 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/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/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 ee8e683..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 - ---------- - METRIC : string - The type of similarity metric being used. - DESCRIPTION : string - A meaningful description of the model used, with references where - appropriate. - """ - - MODEL=None - DESCRIPTION=None - - - 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 4f248f0..0000000 Binary files a/register/samplers/build/lib.macosx-10.6-x86_64-2.7/libsampler.so and /dev/null differ 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 0f2d69c..0000000 Binary files a/register/samplers/build/temp.macosx-10.6-x86_64-2.7/libsampler.o and /dev/null differ 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():