diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 970c6b8..f9654d8 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -1,8 +1,8 @@ - Nathan Faggian - Project coordination, methods development. + Australian Bureau of Meteorology - Stefan Van Der Walt - Project coordination, methods development. + UC Berkley - Riaan Van Den Dool - Methods development \ 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/bento.info b/bento.info new file mode 100644 index 0000000..ce083ef --- /dev/null +++ b/bento.info @@ -0,0 +1,35 @@ +Name: imreg +Version: 0.1 +Summary: Image registration algorithms +DownloadUrl: http://github.com/pyimreg/pyimreg +Description: Image registration algorithms + + Affine, projective, spline, etc. + +Maintainer: Nathan Faggian +MaintainerEmail: N.Faggian@bom.gov.au +License: Apache 2 +Classifiers: + Development Status :: 4 - Beta, + Environment :: Console, + Intended Audience :: Developers, + Intended Audience :: Science/Research, + License :: OSI Approved :: BSD License, + Programming Language :: C, + Programming Language :: Python, + Programming Language :: Python :: 3, + Topic :: Scientific/Engineering, + Operating System :: Microsoft :: Windows, + Operating System :: POSIX, + Operating System :: Unix, + Operating System :: MacOS + +HookFile: bscript +UseBackends: Waf + +Library: + Packages: + imreg, imreg.samplers + Extension: imreg.samplers.libsampler + Sources: + imreg/samplers/libsampler.cpp diff --git a/bscript b/bscript new file mode 100644 index 0000000..2cccc18 --- /dev/null +++ b/bscript @@ -0,0 +1,19 @@ +import sys + +from bento.commands import hooks + +from numpy.distutils.misc_util \ + import \ + get_numpy_include_dirs + + +@hooks.post_configure +def post_configure(context): + conf = context.waf_context + conf.env.INCLUDES = get_numpy_include_dirs() + + if sys.platform == "win32": + conf.options.check_cxx_compiler = "msvc" + else: + conf.options.check_cxx_compiler = "g++" + conf.load("compiler_cxx") diff --git a/documentation/Makefile b/documentation/Makefile new file mode 100644 index 0000000..4156f36 --- /dev/null +++ b/documentation/Makefile @@ -0,0 +1,130 @@ +# Makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +PAPER = +BUILDDIR = _build + +# Internal variables. +PAPEROPT_a4 = -D latex_paper_size=a4 +PAPEROPT_letter = -D latex_paper_size=letter +ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . + +.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest + +help: + @echo "Please use \`make ' where is one of" + @echo " html to make standalone HTML files" + @echo " dirhtml to make HTML files named index.html in directories" + @echo " singlehtml to make a single large HTML file" + @echo " pickle to make pickle files" + @echo " json to make JSON files" + @echo " htmlhelp to make HTML files and a HTML help project" + @echo " qthelp to make HTML files and a qthelp project" + @echo " devhelp to make HTML files and a Devhelp project" + @echo " epub to make an epub" + @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter" + @echo " latexpdf to make LaTeX files and run them through pdflatex" + @echo " text to make text files" + @echo " man to make manual pages" + @echo " changes to make an overview of all changed/added/deprecated items" + @echo " linkcheck to check all external links for integrity" + @echo " doctest to run all doctests embedded in the documentation (if enabled)" + +clean: + -rm -rf $(BUILDDIR)/* + +html: + $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html + @echo + @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." + +dirhtml: + $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml + @echo + @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml." + +singlehtml: + $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml + @echo + @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml." + +pickle: + $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle + @echo + @echo "Build finished; now you can process the pickle files." + +json: + $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json + @echo + @echo "Build finished; now you can process the JSON files." + +htmlhelp: + $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp + @echo + @echo "Build finished; now you can run HTML Help Workshop with the" \ + ".hhp project file in $(BUILDDIR)/htmlhelp." + +qthelp: + $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp + @echo + @echo "Build finished; now you can run "qcollectiongenerator" with the" \ + ".qhcp project file in $(BUILDDIR)/qthelp, like this:" + @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/python-register.qhcp" + @echo "To view the help file:" + @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/python-register.qhc" + +devhelp: + $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp + @echo + @echo "Build finished." + @echo "To view the help file:" + @echo "# mkdir -p $$HOME/.local/share/devhelp/python-register" + @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/python-register" + @echo "# devhelp" + +epub: + $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub + @echo + @echo "Build finished. The epub file is in $(BUILDDIR)/epub." + +latex: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo + @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." + @echo "Run \`make' in that directory to run these through (pdf)latex" \ + "(use \`make latexpdf' here to do that automatically)." + +latexpdf: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo "Running LaTeX files through pdflatex..." + make -C $(BUILDDIR)/latex all-pdf + @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." + +text: + $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text + @echo + @echo "Build finished. The text files are in $(BUILDDIR)/text." + +man: + $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man + @echo + @echo "Build finished. The manual pages are in $(BUILDDIR)/man." + +changes: + $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes + @echo + @echo "The overview file is in $(BUILDDIR)/changes." + +linkcheck: + $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck + @echo + @echo "Link check complete; look for any errors in the above output " \ + "or in $(BUILDDIR)/linkcheck/output.txt." + +doctest: + $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest + @echo "Testing of doctests in the sources finished, look at the " \ + "results in $(BUILDDIR)/doctest/output.txt." diff --git a/documentation/conf.py b/documentation/conf.py new file mode 100644 index 0000000..cf31a96 --- /dev/null +++ b/documentation/conf.py @@ -0,0 +1,216 @@ +# -*- coding: utf-8 -*- +# +# python-register documentation build configuration file, created by +# sphinx-quickstart on Sun Aug 21 10:47:21 2011. +# +# This file is execfile()d with the current directory set to its containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +import sys, os + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +#sys.path.insert(0, os.path.abspath('.')) + +# -- General configuration ----------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +#needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be extensions +# coming with Sphinx (named 'sphinx.ext.*') or your custom ones. +extensions = ['sphinx.ext.autodoc', 'sphinx.ext.autosummary', 'numpydoc'] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix of source filenames. +source_suffix = '.rst' + +# The encoding of source files. +#source_encoding = 'utf-8-sig' + +# The master toctree document. +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' + +# The version info for the project you're documenting, acts as replacement for +# |version| and |release|, also used in various other places throughout the +# built documents. +# +# The short X.Y version. +version = '0.0' +# The full version, including alpha/beta/rc tags. +release = '0.0' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +#language = None + +# There are two options for replacing |today|: either, you set today to some +# non-false value, then it is used: +#today = '' +# Else, today_fmt is used as the format for a strftime call. +#today_fmt = '%B %d, %Y' + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +exclude_patterns = ['_build'] + +# The reST default role (used for this markup: `text`) to use for all documents. +#default_role = None + +# If true, '()' will be appended to :func: etc. cross-reference text. +#add_function_parentheses = True + +# If true, the current module name will be prepended to all description +# unit titles (such as .. function::). +#add_module_names = True + +# If true, sectionauthor and moduleauthor directives will be shown in the +# output. They are ignored by default. +#show_authors = False + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# A list of ignored prefixes for module index sorting. +#modindex_common_prefix = [] + + +# -- Options for HTML output --------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +html_theme = 'default' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +#html_theme_options = {} + +# Add any paths that contain custom themes here, relative to this directory. +#html_theme_path = [] + +# The name for this set of Sphinx documents. If None, it defaults to +# " v documentation". +#html_title = None + +# A shorter title for the navigation bar. Default is the same as html_title. +#html_short_title = None + +# The name of an image file (relative to this directory) to place at the top +# of the sidebar. +#html_logo = None + +# The name of an image file (within the static path) to use as favicon of the +# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 +# pixels large. +#html_favicon = None + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, +# using the given strftime format. +#html_last_updated_fmt = '%b %d, %Y' + +# If true, SmartyPants will be used to convert quotes and dashes to +# typographically correct entities. +#html_use_smartypants = True + +# Custom sidebar templates, maps document names to template names. +#html_sidebars = {} + +# Additional templates that should be rendered to pages, maps page names to +# template names. +#html_additional_pages = {} + +# If false, no module index is generated. +#html_domain_indices = True + +# If false, no index is generated. +#html_use_index = True + +# If true, the index is split into individual pages for each letter. +#html_split_index = False + +# If true, links to the reST sources are added to the pages. +#html_show_sourcelink = True + +# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. +#html_show_sphinx = True + +# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. +#html_show_copyright = True + +# If true, an OpenSearch description file will be output, and all pages will +# contain a tag referring to it. The value of this option must be the +# base URL from which the finished HTML is served. +#html_use_opensearch = '' + +# This is the file name suffix for HTML files (e.g. ".xhtml"). +#html_file_suffix = None + +# Output file base name for HTML help builder. +htmlhelp_basename = 'python-registerdoc' + + +# -- Options for LaTeX output -------------------------------------------------- + +# The paper size ('letter' or 'a4'). +#latex_paper_size = 'letter' + +# The font size ('10pt', '11pt' or '12pt'). +#latex_font_size = '10pt' + +# 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', + u'Nathan Faggian, Stefan Van Der Walt, Riaan Van Den Dool', 'manual'), +] + +# The name of an image file (relative to this directory) to place at the top of +# the title page. +#latex_logo = None + +# For "manual" documents, if this is true, then toplevel headings are parts, +# not chapters. +#latex_use_parts = False + +# If true, show page references after internal links. +#latex_show_pagerefs = False + +# If true, show URL addresses after external links. +#latex_show_urls = False + +# Additional stuff for the LaTeX preamble. +#latex_preamble = '' + +# Documents to append as an appendix to all manuals. +#latex_appendices = [] + +# If false, no module index is generated. +#latex_domain_indices = True + + +# -- Options for manual page output -------------------------------------------- + +# 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', + [u'Nathan Faggian, Stefan Van Der Walt, Riaan Van Den Dool'], 1) +] diff --git a/documentation/deformation_models.rst b/documentation/deformation_models.rst new file mode 100644 index 0000000..a778e85 --- /dev/null +++ b/documentation/deformation_models.rst @@ -0,0 +1,6 @@ +============================= +Deformation models +============================= + +.. automodule:: register.models.model + :members: diff --git a/documentation/index.rst b/documentation/index.rst new file mode 100644 index 0000000..2b85489 --- /dev/null +++ b/documentation/index.rst @@ -0,0 +1,32 @@ +.. python-register documentation master file + +python-register +================ + +Linear and non-linear image registration methods, following the methods defined by: + +| Lucas-Kanade 20 Years On: A Unifying Framework +| Simon Baker, Iain Matthews (2004) +| International Journal of Computer Vision 56 (3) p. 221-255 + +| Fast parametric elastic image registration. +| Jan Kybic, Michael Unser (2003) +| IEEE Transactions on Image Processing 12 (11) p. 1427-1442 + + +Contents: + +.. toctree:: + :maxdepth: 10 + + Deformation models + Similarity metrics + Sampling functions + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` + diff --git a/documentation/make.bat b/documentation/make.bat new file mode 100644 index 0000000..6ee046f --- /dev/null +++ b/documentation/make.bat @@ -0,0 +1,170 @@ +@ECHO OFF + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set BUILDDIR=_build +set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% . +if NOT "%PAPER%" == "" ( + set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS% +) + +if "%1" == "" goto help + +if "%1" == "help" ( + :help + echo.Please use `make ^` where ^ is one of + echo. html to make standalone HTML files + echo. dirhtml to make HTML files named index.html in directories + echo. singlehtml to make a single large HTML file + echo. pickle to make pickle files + echo. json to make JSON files + echo. htmlhelp to make HTML files and a HTML help project + echo. qthelp to make HTML files and a qthelp project + echo. devhelp to make HTML files and a Devhelp project + echo. epub to make an epub + echo. latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter + echo. text to make text files + echo. man to make manual pages + echo. changes to make an overview over all changed/added/deprecated items + echo. linkcheck to check all external links for integrity + echo. doctest to run all doctests embedded in the documentation if enabled + goto end +) + +if "%1" == "clean" ( + for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i + del /q /s %BUILDDIR%\* + goto end +) + +if "%1" == "html" ( + %SPHINXBUILD% -b html %ALLSPHINXOPTS% %BUILDDIR%/html + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The HTML pages are in %BUILDDIR%/html. + goto end +) + +if "%1" == "dirhtml" ( + %SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% %BUILDDIR%/dirhtml + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The HTML pages are in %BUILDDIR%/dirhtml. + goto end +) + +if "%1" == "singlehtml" ( + %SPHINXBUILD% -b singlehtml %ALLSPHINXOPTS% %BUILDDIR%/singlehtml + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The HTML pages are in %BUILDDIR%/singlehtml. + goto end +) + +if "%1" == "pickle" ( + %SPHINXBUILD% -b pickle %ALLSPHINXOPTS% %BUILDDIR%/pickle + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can process the pickle files. + goto end +) + +if "%1" == "json" ( + %SPHINXBUILD% -b json %ALLSPHINXOPTS% %BUILDDIR%/json + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can process the JSON files. + goto end +) + +if "%1" == "htmlhelp" ( + %SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% %BUILDDIR%/htmlhelp + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can run HTML Help Workshop with the ^ +.hhp project file in %BUILDDIR%/htmlhelp. + goto end +) + +if "%1" == "qthelp" ( + %SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% %BUILDDIR%/qthelp + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can run "qcollectiongenerator" with the ^ +.qhcp project file in %BUILDDIR%/qthelp, like this: + echo.^> qcollectiongenerator %BUILDDIR%\qthelp\python-register.qhcp + echo.To view the help file: + echo.^> assistant -collectionFile %BUILDDIR%\qthelp\python-register.ghc + goto end +) + +if "%1" == "devhelp" ( + %SPHINXBUILD% -b devhelp %ALLSPHINXOPTS% %BUILDDIR%/devhelp + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. + goto end +) + +if "%1" == "epub" ( + %SPHINXBUILD% -b epub %ALLSPHINXOPTS% %BUILDDIR%/epub + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The epub file is in %BUILDDIR%/epub. + goto end +) + +if "%1" == "latex" ( + %SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; the LaTeX files are in %BUILDDIR%/latex. + goto end +) + +if "%1" == "text" ( + %SPHINXBUILD% -b text %ALLSPHINXOPTS% %BUILDDIR%/text + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The text files are in %BUILDDIR%/text. + goto end +) + +if "%1" == "man" ( + %SPHINXBUILD% -b man %ALLSPHINXOPTS% %BUILDDIR%/man + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The manual pages are in %BUILDDIR%/man. + goto end +) + +if "%1" == "changes" ( + %SPHINXBUILD% -b changes %ALLSPHINXOPTS% %BUILDDIR%/changes + if errorlevel 1 exit /b 1 + echo. + echo.The overview file is in %BUILDDIR%/changes. + goto end +) + +if "%1" == "linkcheck" ( + %SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% %BUILDDIR%/linkcheck + if errorlevel 1 exit /b 1 + echo. + echo.Link check complete; look for any errors in the above output ^ +or in %BUILDDIR%/linkcheck/output.txt. + goto end +) + +if "%1" == "doctest" ( + %SPHINXBUILD% -b doctest %ALLSPHINXOPTS% %BUILDDIR%/doctest + if errorlevel 1 exit /b 1 + echo. + echo.Testing of doctests in the sources finished, look at the ^ +results in %BUILDDIR%/doctest/output.txt. + goto end +) + +:end diff --git a/documentation/sampling_functions.rst b/documentation/sampling_functions.rst new file mode 100644 index 0000000..4ae0289 --- /dev/null +++ b/documentation/sampling_functions.rst @@ -0,0 +1,9 @@ +======================== +Image resampling methods +======================== + +.. automodule:: register.samplers.sampler + :members: + +.. automodule:: register.samplers.libsampler + :members: diff --git a/documentation/similarity_metrics.rst b/documentation/similarity_metrics.rst new file mode 100644 index 0000000..890129c --- /dev/null +++ b/documentation/similarity_metrics.rst @@ -0,0 +1,6 @@ +======================== +Similarity metrics +======================== + +.. automodule:: register.metrics.metric + :members: diff --git a/examples/README.rst b/examples/README.rst index 1909509..cfeaa66 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -19,17 +19,3 @@ Demonstrates uses of the sci-kit: F : floating image ( non-deformed image) p : warp parameters - + nonlinreg - - Uses an cubic spline deformation model to defrom the image and then (attempts to) minimize: - - || T - W(F;p) || + a*||p|| - - where: - - T : target image (lena deformed) - F : floating image ( non-deformed image) - p : warp parameters - a : regularization term - determines smoothness of the warping. - - diff --git a/examples/haar_features.py b/examples/haar_features.py deleted file mode 100644 index f4ef46f..0000000 --- a/examples/haar_features.py +++ /dev/null @@ -1,27 +0,0 @@ -""" -Detects haar salient features in an image - - -""" - -import scipy.ndimage as nd -from matplotlib.pyplot import imread, plot, imshow, show - -from register.features.detector import detect, HaarDetector - -# Load the image. -image = imread('data/cameraman.png') -#image = nd.zoom(image, 0.50) - -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 - -features = detect(image, HaarDetector, options) - -imshow(image, cmap='gray') - -for id, point in features['points'].items(): - plot(point[1], point[0], 'or') - -show() diff --git a/examples/linreg.py b/examples/linreg.py index 9cb78f3..b1a5732 100644 --- a/examples/linreg.py +++ b/examples/linreg.py @@ -1,4 +1,4 @@ -""" +""" Estimates a linear warp field, the target is a transformed version of lenna: http://en.wikipedia.org/wiki/Lenna @@ -7,16 +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 import register, model, metric +from imreg.samplers import sampler # Form some test data (lena, lena rotated 20 degrees) image = misc.lena() -image = nd.zoom(image, 0.20) template = nd.rotate(image, 20, reshape=False) # Form the affine registration instance. @@ -27,19 +22,12 @@ ) # Coerce the image data into RegisterData. -image = register.RegisterData(image) -template = register.RegisterData(template) - -# Smooth the template and image. -image.smooth(1.5) -template.smooth(1.5) +image = register.RegisterData(image).downsample(2) +template = register.RegisterData(template).downsample(2) # Register. -p, warp, img, error = affine.register( +step, search = affine.register( image, template, - plotCB=plot.gridPlot, - verbose=True + verbose=True, ) - -plot.show() diff --git a/examples/linreg_pyramid.py b/examples/linreg_pyramid.py new file mode 100644 index 0000000..3463b91 --- /dev/null +++ b/examples/linreg_pyramid.py @@ -0,0 +1,49 @@ +""" +Estimates a linear warp field, using an image pyramid to speed up the +computation. +""" + +import scipy.ndimage as nd +import scipy.misc as misc + +from imreg import model, metric, register +from imreg.samplers import sampler + + +# Form some test data (lena, lena rotated 20 degrees) +image = register.RegisterData(misc.lena()) +template = register.RegisterData( + nd.rotate(image.data, 20, reshape=False) + ) + +# Form the registrator. +affine = register.Register( + model.Affine, + metric.Residual, + sampler.CubicConvolution + ) + +fullSearch = [] + +# Image pyramid registration can be executed like so: +pHat = None +for factor in [30., 20. , 10., 5., 2., 1.]: + + if pHat is not None: + scale = downImage.coords.spacing / factor + # FIXME: Find a nicer way to do this. + pHat = model.Affine.scale(pHat, scale) + + downImage = image.downsample(factor) + downTemplate = template.downsample(factor) + + step, search = affine.register( + downImage, + downTemplate, + p=pHat, + verbose=True + ) + + pHat = step.p + + fullSearch.extend(search) diff --git a/examples/nonlinreg.py b/examples/nonlinreg.py deleted file mode 100644 index 4ac518e..0000000 --- a/examples/nonlinreg.py +++ /dev/null @@ -1,77 +0,0 @@ -""" -Estimates a non-linear warp field, the lenna image is randomly deformed using -the spline deformation model. -""" - -import numpy as np -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 - -def warp(image): - """ - Randomly warps an input image using a cubic spline deformation. - """ - coords = register.Coordinates( - [0, image.shape[0], 0, image.shape[1]] - ) - - spline_model = model.Spline(coords) - spline_sampler = sampler.Spline(coords) - - p = spline_model.identity - #TODO: Understand the effect of parameter magnitude: - p += np.random.rand(p.shape[0]) * 100 - 50 - - return spline_sampler.f(image, spline_model.warp(p)).reshape(image.shape) - - -image = misc.lena() -image = nd.zoom(image, 0.30) -template = warp(image) - -# Coerce the image data into RegisterData. -image = register.RegisterData(image) -template = register.RegisterData(template) - -# Smooth the template and image. -image.smooth(0.5) -template.smooth(0.5) - -# Form the affine registration instance. -affine = register.Register( - model.Affine, - metric.Residual, - sampler.CubicConvolution - ) -# Form the spline registration instance. -spline = register.Register( - model.Spline, - metric.Residual, - sampler.CubicConvolution - ) - -# Compute an affine registration between the template and image. -p, warp, img, error = affine.register( - image, - template, - plotCB=plot.gridPlot - ) - -# Compute a nonlinear (spline) registration, initialized with the warp field -# found using the affine registration. -p, warp, img, error = spline.register( - image, - template, - warp=warp, - verbose=True, - plotCB=plot.gridPlot - ) - -plot.show() diff --git a/examples/nonlinreg_image.py b/examples/nonlinreg_image.py deleted file mode 100644 index 1128c5a..0000000 --- a/examples/nonlinreg_image.py +++ /dev/null @@ -1,48 +0,0 @@ -""" -Estimates a linear warp field, between two images - - -The target is a "smile" the image is a frown and the goal is to estimate -the warp field between them. - -The deformation is not manufactured by the spline model, and is a good -(realistic) test of the spline deformation model. - -""" - -from matplotlib.pyplot import imread - -from register.models import model -from register.metrics import metric -from register.samplers import sampler - -from register.visualize import plot -from register import register - -# Form some test data (lena, lena rotated 20 degrees) -image = imread('data/frown.png')[:, :, 0] -template = imread('data/smile.png')[:, :, 0] - -# Form the affine registration instance. -affine = register.Register( - model.Spline, - metric.Residual, - sampler.CubicConvolution - ) - -# Coerce the image data into RegisterData. -image = register.RegisterData(image) -template = register.RegisterData(template) - -# Smooth the template and image. -image.smooth(1.5) -template.smooth(1.5) - -# Register. -p, warp, img, error = affine.register( - image, - template, - plotCB=plot.gridPlot, - verbose=True - ) - -plot.show() diff --git a/examples/nonlinreg_image_features.py b/examples/nonlinreg_image_features.py index bdf68e2..a7363ab 100644 --- a/examples/nonlinreg_image_features.py +++ b/examples/nonlinreg_image_features.py @@ -1,24 +1,22 @@ -""" -Estimates a linear warp field, between two images - +""" +Estimates a warp field, between two images using only image features. The target is a "smile" the image is a frown and the goal is to estimate the warp field between them. - -The deformation is estimated using thin plate splines. - + +The deformation is estimated using thin plate splines and an example of how to +define a custom kernel is shown. """ import numpy as np -import yaml +import yaml import matplotlib.pyplot as plt - -from register.samplers import sampler -from register.visualize import plot -from register import register +from imreg import model, register +from imreg.samplers import sampler -# Load the image and feature data. +# Load the image and feature data. image = register.RegisterData( np.average(plt.imread('data/frown.png'), axis=2), features=yaml.load(open('data/frown.yaml')) @@ -28,23 +26,48 @@ features=yaml.load(open('data/smile.yaml')) ) -# Define a gaussian kernel. -def gaussKernel(r): - var = 50 - return np.exp( -np.power(r,2)/(2*var**2) ) +############################################################################### +# Using the implementation of thin-plate splines. +############################################################################### + +# Form the feature registrator. +feature = register.FeatureRegister( + model=model.ThinPlateSpline, + sampler=sampler.Spline, + ) + +# Perform the registration. +p, warp, img, error = feature.register( + image, + template + ) + +print "Thin-plate Spline kernel error: {}".format(error) + +############################################################################### +# Defining a custom model and registering features. +############################################################################### + +class GaussSpline(model.ThinPlateSpline): + def __init__(self, coordinates): + model.ThinPlateSpline.__init__(self, coordinates) -# Form the tps registration instance. -spline = register.SplineRegister( + def U(self, r): + # Define a gaussian kernel. + var = 5.0 + return np.exp( -np.power(r,2)/(2*var**2) ) + + +# Form feature registrator. +feature = register.FeatureRegister( + model=GaussSpline, sampler=sampler.Spline, - kernel=gaussKernel ) - -# Register using features. -warp, img = spline.register( + +# Perform the registration. +p, warp, img, error = feature.register( image, - template, - vectorized=True, + template ) -plot.featurePlot(image, template, warp, img) -plot.show() \ No newline at end of file +print "Gaussian kernel error: {}".format(error) diff --git a/examples/projective.py b/examples/projective.py new file mode 100644 index 0000000..d2bc7b5 --- /dev/null +++ b/examples/projective.py @@ -0,0 +1,57 @@ +""" +Estimates a projective deformation field, the lenna image is randomly deformed +using the projective deformation model. +""" + +import numpy as np +import scipy.ndimage as nd +import scipy.misc as misc + +from imreg import model, metric, register +from imreg.samplers import sampler + + +def warp(image): + """ + Randomly warps an input image using a projective deformation. + """ + coords = register.Coordinates( + [0, image.shape[0], 0, image.shape[1]] + ) + + mymodel = model.Projective(coords) + mysampler = sampler.CubicConvolution(coords) + + p = mymodel.identity + #TODO: Understand the effect of parameter magnitude: + p += np.random.rand(p.shape[0]) * 0.1 - 0.05 + + return mysampler.f(image, mymodel.warp(p)).reshape(image.shape) + + +image = misc.lena().astype(np.double) +image = nd.zoom(image, 0.30) +template = warp(image) + +# Coerce the image data into RegisterData. +image = register.RegisterData(image) +template = register.RegisterData(template) + +# Smooth the template and image. +image.smooth(0.5) +template.smooth(0.5) + +# Form the affine registration instance. +reg = register.Register( + model.Projective, + metric.Residual, + sampler.CubicConvolution + ) + +# Compute an affine registration between the template and image. +step, search = reg.register( + image, + template, + alpha=0.00002, + verbose=True, + ) diff --git a/register/__init__.py b/imreg/__init__.py similarity index 100% rename from register/__init__.py rename to imreg/__init__.py diff --git a/imreg/metric.py b/imreg/metric.py new file mode 100644 index 0000000..b1fd29b --- /dev/null +++ b/imreg/metric.py @@ -0,0 +1,135 @@ +""" 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/imreg/model.py b/imreg/model.py new file mode 100644 index 0000000..6c20735 --- /dev/null +++ b/imreg/model.py @@ -0,0 +1,770 @@ +""" 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('') + + @property + def identity(self): + raise NotImplementedError('') diff --git a/imreg/register.py b/imreg/register.py new file mode 100644 index 0000000..9485001 --- /dev/null +++ b/imreg/register.py @@ -0,0 +1,531 @@ +""" 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 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/features/__init__.py b/imreg/samplers/__init__.py similarity index 100% rename from register/features/__init__.py rename to imreg/samplers/__init__.py diff --git a/imreg/samplers/libsampler.cpp b/imreg/samplers/libsampler.cpp new file mode 100644 index 0000000..fec46c6 --- /dev/null +++ b/imreg/samplers/libsampler.cpp @@ -0,0 +1,285 @@ +#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/imreg/samplers/ndarray.h similarity index 100% rename from register/samplers/ndarray.h rename to imreg/samplers/ndarray.h diff --git a/register/samplers/numpyctypes.py b/imreg/samplers/numpyctypes.py similarity index 100% rename from register/samplers/numpyctypes.py rename to imreg/samplers/numpyctypes.py diff --git a/register/samplers/sampler.py b/imreg/samplers/sampler.py similarity index 54% rename from register/samplers/sampler.py rename to imreg/samplers/sampler.py index ea11f13..4068269 100644 --- a/register/samplers/sampler.py +++ b/imreg/samplers/sampler.py @@ -5,16 +5,26 @@ 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. - - @param METHOD: the method implemented by the sampler. - @param DESCRIPTION: a meaningful description of the technique used, with - references where appropriate. + 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 @@ -28,11 +38,20 @@ def f(self, array, warp): """ A sampling function, responsible for returning a sampled set of values from the given array. - - @param array: an n-dimensional array (representing an image or volume). - @param coords: array coordinates in cartesian form (n by p). + + 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.') @@ -65,7 +84,7 @@ class Nearest(Sampler): 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) @@ -75,24 +94,90 @@ def f(self, array, warp): """ A sampling function, responsible for returning a sampled set of values from the given array. - - @param array: an n-dimensional array (representing an image or volume). - @param coords: array coordinates in cartesian form (n by p). + + 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) + 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) + 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)' @@ -100,7 +185,7 @@ class CubicConvolution(Sampler): DESCRIPTION=""" Given a coordinate in the array cubic convolution interpolates between 16 (4x4) nearest values. - """ + """ def __init__(self, coordinates): Sampler.__init__(self, coordinates) @@ -110,11 +195,20 @@ def f(self, array, warp): """ A sampling function, responsible for returning a sampled set of values from the given array. - - @param array: an n-dimensional array (representing an image or volume). - @param coords: array coordinates in cartesian form (n by p). + + 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.') @@ -138,8 +232,7 @@ class Spline(Sampler): http://docs.scipy.org/doc/scipy/reference/generated/ scipy.ndimage.interpolation.map_coordinates.html - - """ + """ def __init__(self, coordinates): Sampler.__init__(self, coordinates) @@ -148,11 +241,20 @@ def f(self, array, warp): """ A sampling function, responsible for returning a sampled set of values from the given array. - - @param array: an n-dimensional array (representing an image or volume). - @param coords: array coordinates in cartesian form (n by p). + + 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.') diff --git a/register/samplers/setup.py b/imreg/samplers/setup.py similarity index 100% rename from register/samplers/setup.py rename to imreg/samplers/setup.py diff --git a/register/setup.py b/imreg/setup.py similarity index 59% rename from register/setup.py rename to imreg/setup.py index c558394..32fae68 100644 --- a/register/setup.py +++ b/imreg/setup.py @@ -4,13 +4,8 @@ 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 = Configuration('imreg', parent_package, top_path) config.add_subpackage('samplers') - config.add_subpackage('visualize') - config.add_subpackage('features') return config 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 e37e5d7..0000000 --- a/register/features/haar2d.py +++ /dev/null @@ -1,85 +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): - 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): - 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 c6f99e1..0000000 --- a/register/metrics/metric.py +++ /dev/null @@ -1,93 +0,0 @@ -""" A collection of image similarity metrics. """ - -import numpy as np - - -class Metric(object): - """ - Abstract similarity metric. - - @param METRIC: the type of similarity metric being used. - @param DESCRIPTION: a meaningful description of the metric used, with - references where appropriate. - """ - - METRIC=None - DESCRIPTION=None - - def __init__(self): - pass - - def error(self, warpedImage, template): - """ - Computes the metric. - - @param warpedImage: a numpy array, representing the image. - @param template: a numpy arrary, representing the template. - """ - - raise NotImplementedError('') - - def jacobian(self, model, warpedImage): - """ - Computes the jacobian dP/dE - - @param model: the deformation model. - @param warpedImage: the transformed image. - """ - 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): - """ - Follows the derivations shown in: - - Simon Baker and Iain Matthews. 2004. Lucas-Kanade 20 Years On: A - Unifying Framework. Int. J. Comput. Vision 56, 3 (February 2004). - - @param model: the deformation model. - @param warpedModel: a transformed image. - @return: 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() - - 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): - - 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 6e41f0e..0000000 --- a/register/models/model.py +++ /dev/null @@ -1,304 +0,0 @@ -""" A collection of deformation models. """ - -import numpy as np -import scipy.signal as signal - -class Model(object): - """ - Abstract geometry model. - - @param MODEL: the deformation model used. - @param DESCRIPTION: a meaningful description of the model used, with - references where appropriate. - """ - - MODEL=None - DESCRIPTION=None - - def __init__(self, coordinates): - - self.coordinates = coordinates - - def estimate(self, warp): - """" - Estimates the best fit parameters that define a warp field. - - @param warp: a warp field, representing the warped coordinates. - @return: a set of parameters (n-dimensional array). - """ - raise NotImplementedError('') - - def warp(self, parameters): - """ - Computes the warp field given transformed coordinates. - - @param param: array coordinates of model parameters. - @return: a deformation field. - """ - - coords = self.transform(parameters) - - warp = np.zeros_like(self.coordinates.tensor) - - warp[0] = coords[1].reshape(warp[0].shape) - warp[1] = coords[0].reshape(warp[1].shape) - - # Return the difference warp grid. - return warp - - def transform(self, parameters): - """ - A geometric transformation of coordinates. - - @param param: array coordinates of model parameters. - """ - 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. - """ - - def __init__(self, coordinates): - Model.__init__(self, coordinates) - - @property - def identity(self): - return np.zeros(2) - - def transform(self, parameters): - """ - - Applies an shift transformation to image coordinates. - - @param parameters: a array of shift parameters. - @param coords: array coordinates in cartesian form (n by p). - @return: a transformed set of coordinates. - """ - - T = np.eye(3,3) - T[0,2] = -parameters[0] - T[1,2] = -parameters[1] - - return np.dot(T, self.coordinates.homogenous) - - def jacobian(self): - """ - 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). - - 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. - """ - - def __init__(self, coordinates): - Model.__init__(self, coordinates) - - @property - def identity(self): - return np.zeros(6) - - def transform(self, p): - """ - Applies an affine transformation to image coordinates. - - @param parameters: a array of affine parameters. - @param coords: array coordinates in cartesian form (n by p). - @return: a transformed set of coordinates. - """ - - T = np.array([ - [p[0]+1.0, p[2], p[4]], - [p[1], p[3]+1.0, p[5]], - [0, 0, 1] - ]) - - return np.dot(np.linalg.inv(T), self.coordinates.homogenous) - - def jacobian(self): - """ - 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). - - 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 Spline(Model): - - MODEL='Spline (S)' - - DESCRIPTION=""" - Applies a 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=4): - """ - Follows the derivations in: - - Kybic, J. and Unser, M. (2003). Fast parametric elastic image - registration. IEEE Transactions on Image Processing, 12(11), 1427-1442. - - Computes the spline tensor product and stores the products, as basis - vectors. - - @param order: b-spline order. - @param division: 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. - - @param warp: a warp field, representing the warped coordinates. - @return: a set of parameters (n-dimensional array). - """ - - return np.hstack( - ( - np.dot(np.linalg.pinv(self.basis), - (self.coordinates.tensor[0] - warp[0]).flatten()), - np.dot(np.linalg.pinv(self.basis), - (self.coordinates.tensor[1] - warp[1]).flatten()), - ) - ).T - - - def warp(self, parameters): - """ - Computes the (inverse) warp field given transformed coordinates. - - @param param: array coordinates of model parameters. - @return: a deformation field. - """ - - dwarp = self.transform(parameters) - return self.coordinates.tensor - dwarp - - def transform(self, p): - """ - Applies an spline transformation to image coordinates. - - @param parameters: a array of affine parameters. - @param coords: array coordinates in cartesian form (n by p). - @return: a transformed set of coordinates. - """ - - px = np.array(p[0:self.numberOfParameters]) - py = np.array(p[self.numberOfParameters::]) - - shape = self.coordinates.tensor[0].shape - - return np.array( [ np.dot(self.basis, py).reshape(shape), - np.dot(self.basis, px).reshape(shape) - ] - ) - - def jacobian(self): - """ - Follows the derivations shown in: - - Kybic, J., & Unser, M. (2003). Fast parametric elastic image - registration. IEEE Transactions on Image Processing, 12(11), 1427-1442. - - 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 039fa48..0000000 --- a/register/register.py +++ /dev/null @@ -1,448 +0,0 @@ -""" A top level registration module """ - -import collections -import numpy as np -import scipy.ndimage as nd - - -def _smooth(image, variance): - """ - A simple image smoothing method - using a Gaussian kernel. - @param image: the input image, a numpy ndarray object. - @param variance: the width of the smoothing kernel. - @return: the smoothing input image. - """ - return np.real( - np.fft.ifft2( - nd.fourier_gaussian( - np.fft.fft2(image), - variance - ) - ) - ) - - -class Coordinates(object): - """ - A container for grid coordinates. - """ - def __init__(self, domain, spacing=None): - - 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): - """ - A container for registration data. - """ - def __init__(self, data, coords=None, features=None): - - self.data = data - - if not coords: - self.coords = Coordinates( - [0, data.shape[0], 0, data.shape[1]] - ) - else: - self.coords = coords - - # Features are (as a starting point a dictionary) which define - # labelled salient image coordinates (point features). - - self.features = features - - - def smooth(self, variance): - """ - A simple image smoothing method - using a Gaussian kernel. - @param variance: the width of the smoothing kernel. - """ - self.data = _smooth(self.data, variance) - - -class Register(object): - """ - A registration class for estimating the deformation model parameters that - best solve: - - f( W(I;p), T ) - - where: - f : is a similarity metric. - W(x;p): is a deformation model (defined by the parameter set p). - I : is an input image (to be deformed). - T : is a template (which is a deformed version of the input). - - """ - - # The optimization step cache. - optStep = collections.namedtuple('optStep', 'error p deltaP') - - # The maximum number of optimization iterations. - MAX_ITER = 200 - - # The maximum numver of bad (incorrect) optimization steps. - MAX_BAD = 20 - - def __init__(self, model, metric, sampler): - - self.model = model - self.metric = metric - self.sampler = sampler - - def __deltaP(self, J, e, alpha, p=None): - """ - Compute the parameter update. - - Refer to the Levernberg-Marquardt algorithm: - http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm - - @param J: dE/dP the relationship between image differences and model - parameters. - @param e: the difference between the image and template. - @param alpha: the dampening factor. - @keyword p: the current parameter set. - @return: deltaP, the set of model parameter updates. (p x 1). - """ - - 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): - """ - Returns the dampening value. - - Refer to the Levernberg-Marquardt algorithm: - http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm - - @param alpha: a dampening factor. - @param decreasing: a boolean indicating that the error function is - decreasing. - @return: an adjusted dampening factor. - """ - return alpha / 10. if decreasing else alpha * 10. - - def register(self, - image, - template, - p=None, - alpha=None, - warp=None, - plotCB=None, - verbose=False): - """ - Performs an image registration. - @param image: the floating image. - @param template: the target image. - @keyword p: a list of parameters, (first guess). - @keyword alpha: the dampening factor. - @keyword warp: the warp field (first guess). - @keyword plotCB: a debug plotting function. - @keyword verbose: a debug flag for text status updates. - """ - - #TODO: Determine the common coordinate system. - # if image.coords != template.coords: - # 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 warp is not None: - # Estimate p, using the warp field. - p = model.estimate(warp) - - p = model.identity if p is None else p - deltaP = np.zeros_like(p) - - search = [] - alpha = alpha if alpha is not None else 1e-4 - decreasing = True - badSteps = 0 - - for itteration in range(0,self.MAX_ITER): - - # Compute the warp field (warp field is the inverse warp) - warp = model.warp(p) - - # Sample the image using the inverse warp. - warpedImage = _smooth( - sampler.f(image.data, warp).reshape(image.data.shape), - 0.5 - ) - - # Evaluate the error metric. - e = metric.error(warpedImage, template.data) - - searchStep = self.optStep(error=np.abs(e).sum(), - p=p, - deltaP=deltaP, - ) - - if (len(search) > 1): - - decreasing = (searchStep.error < search[-1].error) - - alpha = self.__dampening( - alpha, - decreasing - ) - - if decreasing: - - 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 iteration. - p = search[-1].p - continue - - # Computes the derivative of the error with respect to model - # parameters. - - J = metric.jacobian(model, warpedImage) - - deltaP = self.__deltaP( - J, - e, - alpha, - p=p - ) - - # Evaluate stopping condition: - if np.dot(deltaP.T, deltaP) < 1e-4: - break - - p += deltaP - - if verbose and decreasing: - print ('{0}\n' - 'iteration : {1} \n' - 'parameters : {2} \n' - 'error : {3} \n' - '{0}\n' - ).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) - - return p, warp, warpedImage, searchStep.error - - -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): - """ - Compute the parameter update. - """ - - 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): - """ - Returns the dampening value, without adjustment. - - @param alpha: a dampening factor. - @param decreasing: a boolean indicating that the error function is - decreasing. - @return: an adjusted dampening factor. - """ - return alpha - - -class SplineRegister(): - """ - A registration class for estimating the deformation field which minimizes - feature differences using a thin-plate-spline interpolant. - """ - - def __init__(self, sampler, kernel=None): - - self.sampler = sampler - self.kernel = kernel if kernel is not None else None - - def U(self, r): - """ - This is a kernel function applied to solve the biharmonic equation. - @param r: - """ - - if not self.kernel: - return np.multiply( -np.power(r,2), np.log(np.power(r,2) + 1e-20)) - else: - return self.kernel(r) - - ## Gaussian kernel - ##var = 5.0 - ##return np.exp( -pow(r,2)/(2*var**2) ) - - def approximate(self, p0, p1): - """ - Approximates the thinplate spline coefficients, following derivations - shown 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. - """ - - 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))) ) - - Y = np.matrix(Y) - - Linv = np.matrix(np.linalg.inv(L)) - - return L, Linv*Y - - def register(self, - image, - template, - vectorized=True): - """ - Performs an image (feature based) registration. - - @param image: a floating image, registerData object. - @param template: a target image, registerData object. - """ - - 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 image and template features to register.') - - # Note the inverse warp is estimated here. - - p0 = np.array(templatePoints) - p1 = np.array(imagePoints) - - _L, model = self.approximate(p0, p1) - - # For all coordinates in the register data, evaluate the - # thin-plate-spline. - - warp = np.zeros_like(image.coords.tensor) - - affine = model[-3:, :] - weights = model[:-3, :] - - if vectorized: - - # Vectorized extrapolation, looping through arrays in python is - # slow therefore wherever possible attempt to unroll loops. Below - # is an example: - - Xvec = np.matrix(image.coords.tensor[1].flatten(0)).T - Yvec = np.matrix(image.coords.tensor[0].flatten(0)).T - - # Fast matrix multiplication approach: - Px = np.matrix(np.tile(p0[:,0], (Xvec.shape[0], 1))) - Wx = np.matrix(np.tile(weights[:,0].T, (Xvec.shape[0], 1))) - Bx = np.matrix(np.tile(Xvec, (1, p0.shape[0]))) - Ax = np.matrix(np.tile(affine[:,0].T, (Xvec.shape[0], 1))) - - Py = np.matrix(np.tile(p0[:,1], (Xvec.shape[0], 1))) - Wy = np.matrix(np.tile(weights[:,1].T, (Xvec.shape[0], 1))) - By = np.matrix(np.tile(Yvec, (1, p0.shape[0]))) - Ay = np.matrix(np.tile(affine[:,1].T, (Xvec.shape[0], 1))) - - # Form the R matrix: - R = self.U( np.sqrt( np.power(Px - Bx,2) + np.power(Py - By,2)) ) - - # Compute the sum of the weighted R matrix, row wise. - Rx = np.sum( np.multiply(Wx, R), 1 ) - Ry = np.sum( np.multiply(Wy, R), 1 ) - - one = np.ones_like(Xvec) - a = np.hstack(( one, Xvec, Yvec, one)) - - warp[1] = np.sum( np.multiply(a, np.hstack((Ax, Rx)) ), 1 ).reshape(warp[1].shape) - warp[0] = np.sum( np.multiply(a, np.hstack((Ay, Ry)) ), 1 ).reshape(warp[0].shape) - - else: - - # Slow nested loop approach: - for x in xrange(0, image.coords.domain[1]): - for y in xrange(0, image.coords.domain[3]): - - # Refer to page 570 (of BookStein paper) first column, last - # equation (relating map coordinates to coefficients). - - zx = 0.0 - zy = 0.0 - - for n in range(0, len(p0[:,0])): - r = np.sqrt( (p0[n,0] - x)**2 + (p0[n,1] - y)**2 ) - zx += float(weights[n,0])*float(self.U(r)) - zy += float(weights[n,1])*float(self.U(r)) - - warp[0][y,x] = affine[0,1] + affine[1,1]*x + affine[2,1]*y + zy - warp[1][y,x] = affine[0,0] + affine[1,0]*x + affine[2,0]*y + zx - - img = sampler.f(image.data, warp).reshape(image.data.shape) - - return warp, img 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 132dd23..0000000 --- a/register/samplers/libsampler.cpp +++ /dev/null @@ -1,112 +0,0 @@ -#include -#include "ndarray.h" -#include "math.h" - -using namespace std; - -extern "C" { - -int nearest(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); - - for (int i = 0; i < image.getShape(0); i++) - { - for (int j = 0; j < image.getShape(1); j++) - { - di = (int)warp[0][i][j]; - dj = (int)warp[1][i][j]; - - if ( ( di < rows && di >= 0 ) && - ( dj < cols && dj >= 0 ) ) - { - result[i][j] = image[di][dj]; - } - else - { - result[i][j] = 0.0; - } - } - } - - 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; -} - -} // end extern "C" diff --git a/register/visualize/__init__.py b/register/visualize/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/register/visualize/plot.py b/register/visualize/plot.py deleted file mode 100644 index ce0015c..0000000 --- a/register/visualize/plot.py +++ /dev/null @@ -1,186 +0,0 @@ -''' (Debug utility) Defines a set of plotting callback functions ''' - -import matplotlib.pyplot as plt -import numpy as np - -IMAGE_ORIGIN=None -IMAGE_COLORMAP='gray' -IMAGE_VMIN=None -IMAGE_VMAX=None - -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.ion() -plt.rcParams.update(params); - -def show(): - - plt.ioff() - plt.show() - -def coordPlt(grid, buffer=10, step=10): - """ - 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..6f30c49 100644 --- a/setup.py +++ b/setup.py @@ -1,18 +1,17 @@ #! /usr/bin/env python -descr = """Place long description here. - +DISTNAME = 'imreg' +DESCRIPTION = 'Image registration toolkit' +LONG_DESCRIPTION = """ +"imreg" is an image registration package for python that makes it easy to +automatically align image data. """ - -DISTNAME = 'python-register' -DESCRIPTION = 'Image registration toolbox for SciPy' -LONG_DESCRIPTION = descr -MAINTAINER = 'Nathan Faggian' -MAINTAINER_EMAIL = 'nathan.faggian@gmail.com' -URL = '' -LICENSE = 'Apache License (2.0)' -DOWNLOAD_URL = '' -VERSION = '0.1' +MAINTAINER = 'Nathan Faggian, Riaan Van Den Dool, Stefan Van Der Walt' +MAINTAINER_EMAIL = 'nathan.faggian@gmail.com' +URL = 'pyimreg.github.com' +LICENSE = 'Apache License (2.0)' +DOWNLOAD_URL = '' +VERSION = '0.1' import os import setuptools @@ -34,8 +33,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 @@ -63,7 +61,6 @@ def configuration(parent_package='', top_path=None): install_requires=[], packages=setuptools.find_packages(), include_package_data=True, - zip_safe=False, # the package can run out of an .egg file - + zip_safe=False, cmdclass={'build_py': build_py}, ) diff --git a/tests/test_detectors.py b/tests/test_detectors.py deleted file mode 100644 index 56aae7a..0000000 --- a/tests/test_detectors.py +++ /dev/null @@ -1,27 +0,0 @@ -import numpy as np -import register.samplers.sampler as sampler -import time -import matplotlib.pyplot as plt -import os - -from register.features.detector import detect, HaarDetector - -def test_haardetector(): - """ - Excersize the Haar feature detector. - Asserts that some basic test cases are correct. - - """ - - path = os.path.dirname(__file__) - image = plt.imread('%s/../examples/data/cameraman.png' % path) - - options = {} - options['levels'] = 5 - options['threshold'] = 0.2 - options['locality'] = 5 - - features = detect(image, HaarDetector, options, debug=True) - - assert len(features['points'].items()) > 0 - diff --git a/tests/test_haar2d.py b/tests/test_haar2d.py deleted file mode 100644 index 45382f4..0000000 --- a/tests/test_haar2d.py +++ /dev/null @@ -1,19 +0,0 @@ -import numpy as np -import register.samplers.sampler as sampler -import time - -from register.features.haar2d import haar2d, ihaar2d - -def test_haar2d(): - """ - Asserts that some basic test cases are correct. - - """ - - assert haar2d(np.random.random([5,3]),2,debug=True).shape == (8,4), "Transform data must be padded to compatible shape." - assert haar2d(np.random.random([8,4]),2,debug=True).shape == (8,4), "Transform data must be padded to compatible shape, only if neccersary." - - image = np.random.random([5,3]) - haart = haar2d(image, 3, debug=False) - haari = ihaar2d(haart, 3, debug=False)[:image.shape[0], :image.shape[1]] - assert (image - haari < 0.0001).all(), "Transform must be circular." diff --git a/tests/test_models.py b/tests/test_models.py new file mode 100644 index 0000000..69d3234 --- /dev/null +++ b/tests/test_models.py @@ -0,0 +1,117 @@ +import numpy as np + +from imreg import model, register +from imreg.samplers import sampler + +import scipy.misc as misc +import scipy.ndimage as nd + + +def test_shift(): + """ + Asserts that the feature point alignment error is sufficiently small. + """ + + # Form a dummy coordinate class. + coords = register.Coordinates( + [0, 10, 0, 10] + ) + + # Form corresponding feature sets. + p0 = np.array([0, 0, 0, 1, 1, 0, 1, 1]).reshape(4, 2) + p1 = p0 + 2.0 + + shift = model.Shift(coords) + + _parameters, error = shift.fit(p0, p1) + + print _parameters + + # Assert that the alignment error is small. + assert error <= 1.0, "Unexpected large alignment error : {} grid units".format(error) + + +def test_affine(): + """ + Asserts that the feature point alignment error is sufficiently small. + """ + + # Form a dummy coordinate class. + coords = register.Coordinates( + [0, 10, 0, 10] + ) + + # Form corresponding feature sets. + p0 = np.array([0, 0, 0, 1, 1, 0, 1, 1]).reshape(4,2) + p1 = p0 + 2.0 + + affine = model.Affine(coords) + + _parameters, error = affine.fit(p0, p1) + + # Assert that the alignment error is small. + assert error <= 1.0, "Unexpected large alignment error : {} grid units".format(error) + + +def test_thinPlateSpline(): + """ + Asserts that the feature point alignment error is sufficiently small. + """ + + # Form a dummy coordinate class. + coords = register.Coordinates( + [0, 10, 0, 10] + ) + + # Form corresponding feature sets. + p0 = np.array([0,1, -1,0, 0,-1, 1, 0]).reshape(4,2) + p1 = np.array([0,0.75, -1, 0.25, 0, -1.25, 1, 0.25]).reshape(4,2) + + spline = model.ThinPlateSpline(coords) + + _parameters, error = spline.fit(p0, p1) + + # Assert that the alignment error is small. + assert error < 1.0, "Unexpected large alignment error." + + +def test_thinPlateSplineApproximate(): + """ + Asserts that the computed K, P, L and V matrices are formed correctly. + + Refer to: + + 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. + """ + + # Form a dummy coordinate class. + coords = register.Coordinates( + [0, 100, 0, 100] + ) + + # Form corresponding feature sets. + p0 = np.array([0,1, -1,0, 0,-1, 1, 0]).reshape(4,2) + p1 = np.array([0,0.75, -1, 0.25, 0, -1.25, 1, 0.25]).reshape(4,2) + + spline = model.ThinPlateSpline(coords) + + _parameters, _error, L = spline.fit(p0, p1, lmatrix=True) + + # This expected L matrix is derived from the symmetric example in the + # referenced paper. + expectedL = np.array( + [[ 0. , -1.3863, -5.5452, -1.3863, 1. , 0. , 1. ], + [-1.3863, 0. , -1.3863, -5.5452, 1. , -1. , 0. ], + [-5.5452, -1.3863, 0. , -1.3863, 1. , 0. , -1. ], + [-1.3863, -5.5452, -1.3863, 0. , 1. , 1. , 0. ], + [ 1. , 1. , 1. , 1. , 0. , 0. , 0. ], + [ 0. , -1. , 0. , 1. , 0. , 0. , 0. ], + [ 1. , 0. , -1. , 0. , 0. , 0. , 0. ]] + ) + + assert np.allclose(L, expectedL), \ + "The expected L matrix was not derived." + + diff --git a/tests/test_register.py b/tests/test_register.py index 2156487..e572239 100644 --- a/tests/test_register.py +++ b/tests/test_register.py @@ -3,25 +3,18 @@ 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 +from imreg import model, metric, register +from imreg.samplers import sampler -from register import register def warp(image, p, model, sampler): """ - Warps an image given a deformation model a set of parameters. - - @param image: an numpy ndarray. - @param p: warp parameters. - @param model: a deformation model. - + Warps an image. """ coords = register.Coordinates( [0, image.shape[0], 0, image.shape[1]] ) - + model = model(coords) sampler = sampler(coords) @@ -104,19 +97,16 @@ def test_shift(image, template, p): # Coerce the image data into RegisterData. image = register.RegisterData(image) template = register.RegisterData(template) - - # Smooth - image.smooth(0.5) - template.smooth(0.5) - - _p, _warp, _img, _error = shift.register( + + + step, _search = shift.register( image, template ) - assert np.allclose(p, _p, atol=0.5), \ + assert np.allclose(p, step.p, atol=0.5), \ "Estimated p: {} not equal to p: {}".format( - _p, + step.p, p ) @@ -135,18 +125,14 @@ def test_affine(image, template, p): # Coerce the image data into RegisterData. image = register.RegisterData(image) template = register.RegisterData(template) - - # Smooth - image.smooth(0.5) - template.smooth(0.5) - - _p, _warp, _img, _error = affine.register( + + step, _search = affine.register( image, template ) - assert np.allclose(p, _p, atol=0.5), \ + assert np.allclose(p, step.p, atol=0.5), \ "Estimated p: {} not equal to p: {}".format( - _p, + step.p, p ) diff --git a/tests/test_register_data.py b/tests/test_register_data.py new file mode 100644 index 0000000..fdc723c --- /dev/null +++ b/tests/test_register_data.py @@ -0,0 +1,14 @@ +import scipy.misc as misc + +from imreg import register + +def test_downsample(): + """ + Tests register data down-sampling. + """ + image = register.RegisterData(misc.lena()) + for factor in [1, 2, 4, 6, 8 ,10]: + subSampled = image.downsample(factor) + assert subSampled.data.shape[0] == image.data.shape[0] / factor + assert subSampled.data.shape[1] == image.data.shape[1] / factor + assert subSampled.coords.spacing == factor diff --git a/tests/test_register_features.py b/tests/test_register_features.py index b0fc7cc..c5ffdcd 100644 --- a/tests/test_register_features.py +++ b/tests/test_register_features.py @@ -1,19 +1,19 @@ import numpy as np -import time -import register.samplers.sampler as sampler -from register import register +from imreg import model, register +from imreg.samplers import sampler def test_register(): - """ + """ Top level registration of a simple unit square. """ + img = np.zeros((100,100)) img[25:75, 25:75] = 1 - + image = register.RegisterData( - img, + img, features={ 'points': { @@ -23,143 +23,32 @@ def test_register(): '004': [75, 75], } } - ) - - template = register.RegisterData( - img, - features={ - 'points': - { - '001': [35, 35], - '002': [35, 85], - '003': [85, 35], - '004': [50, 50], - } - } - ) - - # Form a thinplate spline "registrator" - - spline = register.SplineRegister( - sampler.CubicConvolution ) - - warp, img = spline.register(image, template) - - assert not np.allclose(warp, np.zeros_like(warp)), \ - "Estimated warp field is zero." - - -def test_vectorized(): - """ - Asserts that the execution time (on average) of vectorized code is *faster*. - """ - - # Define some dummy data. - img = np.zeros((100,100)) - img[25:75, 25:75] = 1 - - image = register.RegisterData( - img, - features={ - 'points': - { - '001': [25, 25], - '002': [25, 75], - '003': [75, 25], - '004': [75, 75], - } - } - ) - + template = register.RegisterData( - img, + img, features={ 'points': { '001': [35, 35], '002': [35, 85], '003': [85, 35], - '004': [50, 50], + '004': [70, 70], } } - ) - - # Form a thinplate spline "registrator" - - spline = register.SplineRegister( - sampler.CubicConvolution - ) - - times = np.zeros(10) - for i in range(0,10): - t1 = time.time() - _warp, img = spline.register(image, template) - t2 = time.time() - times[i] = (t2-t1)*1000.0 - - print 'Vectorized : {0}x{0} image - {1:0.3f} ms'.format( - 100, - np.average(times) - ) - - vtimes = np.zeros(10) - for i in range(0,10): - t1 = time.time() - _warp, img = spline.register(image, template, vectorized=False) - t2 = time.time() - vtimes[i] = (t2-t1)*1000.0 - - print 'Untouched : {0}x{0} image - {1:0.3f} ms'.format( - 100, - np.average(vtimes) - ) - - assert np.average(times) < np.average(vtimes), \ - "Vectorized code is slower than non-vectorized code. Not good." - - -def test_approximate(): - """ - Asserts that the computed K, P, L and V matrices are formed correctly. - - Refer to: - - 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. - """ - - # Form a dummy coordinate class. - coords = register.Coordinates( - [0, 100, 0, 100] ) - - # Form corresponding feature sets. - p0 = np.array([0,1, -1,0, 0,-1, 1, 0]).reshape(4,2) - p1 = np.array([0,0.75, -1, 0.25, 0, -1.25, 1, 0.25]).reshape(4,2) - - # Form a thinplate spline "registrator" - - spline = register.SplineRegister( - sampler.CubicConvolution(coords) + + # Form feature registrator. + feature = register.FeatureRegister( + model=model.Shift, + sampler=sampler.Spline, ) - - L, LinvY = spline.approximate(p0, p1) - - # This expected L matrix is derived from the symmetric example in the - # referenced paper. - expectedL = np.array( - [[ 0. , -1.3863, -5.5452, -1.3863, 1. , 0. , 1. ], - [-1.3863, 0. , -1.3863, -5.5452, 1. , -1. , 0. ], - [-5.5452, -1.3863, 0. , -1.3863, 1. , 0. , -1. ], - [-1.3863, -5.5452, -1.3863, 0. , 1. , 1. , 0. ], - [ 1. , 1. , 1. , 1. , 0. , 0. , 0. ], - [ 0. , -1. , 0. , 1. , 0. , 0. , 0. ], - [ 1. , 0. , -1. , 0. , 0. , 0. , 0. ]] + + # Perform the registration. + _p, warp, _img, _error = feature.register( + image, + template ) - - assert np.allclose(L, expectedL), \ - "The expected L matrix was not derived." - - \ No newline at end of file + + assert not np.allclose(warp, np.zeros_like(warp)), \ + "Estimated warp field is zero." diff --git a/tests/test_samplers.py b/tests/test_samplers.py index 14439b3..6fdc63a 100644 --- a/tests/test_samplers.py +++ b/tests/test_samplers.py @@ -1,64 +1,79 @@ -import numpy as np -import register.samplers.sampler as sampler import time -from register import register +import numpy as np + +from imreg.samplers import sampler +from imreg import register + def test_sampler(): """ - Asserts that NN < Cubic < Spline, over a range of image resolutions. - + Asserts that NN < Bilinear < Cubic < Spline, over a range of image resolutions. + If (one day) something really amazing happens and the scipy map_coordiantes method is significantly faster we could favour that as a default. - + """ - - for n in range(128, 2024, 128): + + for n in range(128, 1024, 128): coords = register.Coordinates( [0, n, 0, n] ) - + img = np.random.rand(n,n) warp = np.random.rand(2,n,n) - + # nearest neighbour sampler - ctypes nearest = sampler.Nearest(coords) - + ntimes = np.zeros(10) for i in range(0,10): t1 = time.time() nearest.f(img, warp) t2 = time.time() ntimes[i] = (t2-t1)*1000.0 - + print 'Nearest : {0}x{0} image - {1:0.3f} ms'.format(n, np.average(ntimes)) - - # cubic convolution sampler- ctypes + + # cubic convolution sampler - ctypes + bilinear = sampler.Bilinear(coords) + + btimes = np.zeros(10) + for i in range(0,10): + t1 = time.time() + bilinear.f(img, warp) + t2 = time.time() + btimes[i] = (t2-t1)*1000.0 + + print 'Bilinear : {0}x{0} image - {1:0.3f} ms'.format(n, np.average(btimes)) + + # cubic convolution sampler - ctypes cubic = sampler.CubicConvolution(coords) - + ctimes = np.zeros(10) for i in range(0,10): t1 = time.time() cubic.f(img, warp) t2 = time.time() ctimes[i] = (t2-t1)*1000.0 - + print 'Cubic : {0}x{0} image - {1:0.3f} ms'.format(n, np.average(ctimes)) # spline sampler - scipy buffered? ctypes? spline = sampler.Spline(coords) - + stimes = np.zeros(10) for i in range(0,10): t1 = time.time() spline.f(img, warp) t2 = time.time() stimes[i] = (t2-t1)*1000.0 - + print 'Spline : {0}x{0} image - {1:0.3f} ms'.format(n, np.average(stimes)) print '====================================' - + assert np.average(ntimes) < np.average(ctimes) + assert np.average(ntimes) < np.average(btimes) assert np.average(ntimes) < np.average(stimes) + assert np.average(btimes) < np.average(ctimes) assert np.average(ctimes) < np.average(stimes) -