Merge pull request #488 from odlgroup/linearized_deformations_simplified

Linearized deformations based on displacement vector fields.

doc/source/prs.inc

@@ -52,6 +52,7 @@
 .. _PR 480: https://github.com/odlgroup/odl/pull/480
 .. _PR 482: https://github.com/odlgroup/odl/pull/482
 .. _PR 487: https://github.com/odlgroup/odl/pull/487
+.. _PR 488: https://github.com/odlgroup/odl/pull/488
 .. _PR 489: https://github.com/odlgroup/odl/pull/489
 .. _PR 490: https://github.com/odlgroup/odl/pull/490
 .. _PR 491: https://github.com/odlgroup/odl/pull/491

doc/source/release_notes.rst

@@ -9,6 +9,10 @@ Release Notes
 Next release
 ============
 
+New features
+------------
+- First implementation of a deformation model based on linearized deformations using displacement
+  vector fields. (`PR 488`)
 
 ODL 0.3.1 Release Notes (2016-08-15)
 ====================================

examples/deform/linearized_fixed_displacement.py

@@ -0,0 +1,82 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Example using the operator of fixed-displacement linearized deformation.
+
+The linearized deformation operator with fixed displacement ``v`` maps
+a given template ``I`` to the function ``x --> I(x + v(x))``.
+This operator is linear (in ``I``).
+
+This example considers a 2D case, where the displacement field ``v``
+is a Gaussian in each component, with positive sign in the first and
+negative sign in the second component. Note that in the deformed image,
+the value at ``x`` is **taken from** the original image at ``x + v(x)``,
+hence the values are moved by ``-v(x)`` when comparing deformed and
+original templates.
+
+The adjoint is based on an approximation given by the transformation
+``x --> I(x - v(x))``, multiplied by an integral transformation factor
+``exp(-div(v))``.
+"""
+
+import numpy as np
+import odl
+
+
+# --- Create template and displacement field --- #
+
+
+# Template space: discretized functions on the rectangle [-1, 1]^2 with
+# 100 samples per dimension.
+templ_space = odl.uniform_discr([-1, -1], [1, 1], (100, 100))
+
+# The template is a rectangle of size 1.0 x 0.5
+template = odl.phantom.cuboid(templ_space, [-0.5, -0.25], [0.5, 0.25])
+
+# Create a product space for displacement field
+disp_field_space = templ_space.vector_field_space
+
+# Define a displacement field that bends the template a bit towards the
+# upper left. We use a list of 2 functions and discretize it using the
+# disp_field_space.element() method.
+sigma = 0.5
+disp_func = [
+    lambda x: 0.4 * np.exp(-(x[0] ** 2 + x[1] ** 2) / (2 * sigma ** 2)),
+    lambda x: -0.3 * np.exp(-(x[0] ** 2 + x[1] ** 2) / (2 * sigma ** 2))]
+
+disp_field = disp_field_space.element(disp_func)
+
+# Show template and displacement field
+template.show('Template')
+disp_field.show('Displacement field')
+
+
+# --- Apply LinDeformFixedDisp and its adjoint --- #
+
+
+# Initialize the deformation operator with fixed displacement
+deform_op = odl.deform.LinDeformFixedDisp(disp_field)
+
+# Apply the deformation operator to get the deformed template.
+deformed_template = deform_op(template)
+
+# Evaluate the adjoint at the same template.
+adj_result = deform_op.adjoint(template)
+
+# Show results
+deformed_template.show('Deformed template')
+adj_result.show('Adjoint applied to the template')

examples/deform/linearized_fixed_template.py

@@ -0,0 +1,91 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Example using the operator of fixed-template linearized deformation.
+
+The linearized deformation operator with fixed template (image) ``I`` maps
+a given displacement field ``v`` to the function ``x --> I(x + v(x))``.
+
+This example consider a 2D case, where the displacement field ``v``
+is a Gaussian in each component, with positive sign in the first and
+negative sign in the second component. Note that in the deformed image,
+the value at ``x`` is **taken from** the original image at ``x + v(x)``,
+hence the values are moved by ``-v(x)`` when comparing deformed and
+original templates.
+
+The derivative and its adjoint are based on the deformation of the
+gradient of the template, hence the result is expected to be some kind of
+edge image or "edge vector field", respectively.
+"""
+
+import numpy as np
+import odl
+
+
+# --- Create template and displacement field --- #
+
+
+# Template space: discretized functions on the rectangle [-1, 1]^2 with
+# 100 samples per dimension. Usage of 'linear' interpolation ensures that
+# the template gradient is well-defined.
+templ_space = odl.uniform_discr([-1, -1], [1, 1], (100, 100), interp='linear')
+
+# The template is a rectangle of size 1.0 x 0.5
+template = odl.phantom.cuboid(templ_space, [-0.5, -0.25], [0.5, 0.25])
+
+# Create a product space for displacement field
+disp_field_space = templ_space.vector_field_space
+
+# Define a displacement field that bends the template a bit towards the
+# upper left. We use a list of 2 functions and discretize it using the
+# disp_field_space.element() method.
+sigma = 0.5
+disp_func = [
+    lambda x: 0.4 * np.exp(-(x[0] ** 2 + x[1] ** 2) / (2 * sigma ** 2)),
+    lambda x: -0.3 * np.exp(-(x[0] ** 2 + x[1] ** 2) / (2 * sigma ** 2))]
+
+disp_field = disp_field_space.element(disp_func)
+
+# Show template and displacement field
+template.show('Template')
+disp_field.show('Displacement field')
+
+
+# --- Apply LinDeformFixedTempl, derivative and its adjoint --- #
+
+
+# Initialize the deformation operator with fixed template
+deform_op = odl.deform.LinDeformFixedTempl(template)
+
+# Apply the deformation operator to get the deformed template.
+deformed_template = deform_op(disp_field)
+
+# Initialize the derivative of the deformation operator at the
+# given displacement field. The result is again an operator.
+deform_op_deriv = deform_op.derivative(disp_field)
+
+# Evaluate the derivative at the vector field that has value 1 everywhere,
+# i.e. the global shift by (-1, -1).
+deriv_result = deform_op_deriv(disp_field_space.one())
+
+# Evaluate the adjoint of derivative at the image that is 1 everywhere.
+deriv_adj_result = deform_op_deriv.adjoint(templ_space.one())
+
+# Show results
+deformed_template.show('Deformed template')
+deriv_result.show('Operator derivative applied to one()')
+deriv_adj_result.show('Adjoint of the derivative applied to one()')

odl/__init__.py

@@ -27,7 +27,7 @@
 
 __version__ = '0.3.1'
 __all__ = ('diagnostics', 'discr', 'operator', 'set', 'space', 'solvers',
-           'tomo', 'trafos', 'util', 'phantom')
+           'tomo', 'trafos', 'util', 'phantom', 'deform')
 
 # Propagate names defined in __all__ of all submodules into the top-level
 # module
@@ -54,6 +54,7 @@
 from . import tomo
 from . import util
 from . import phantom
+from . import deform
 
 from .util import test
 __all__ += ('test',)

odl/deform/__init__.py

@@ -0,0 +1,27 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Operators and functions for linearized deformations in ODL."""
+
+from __future__ import absolute_import
+
+
+__all__ = ()
+
+
+from .linearized import *
+__all__ += linearized.__all__