Add support for tensorflow together with ODL

odl/__init__.py

@@ -45,6 +45,5 @@
 from . import datasets
 from . import contrib
 
-
 from .util import test
 __all__ += ('test',)

odl/contrib/fom/__init__.py

@@ -12,3 +12,6 @@
 
 from .supervised import *
 __all__ += supervised.__all__
+
+from .unsupervised import *
+__all__ += unsupervised.__all__

odl/contrib/fom/supervised.py

@@ -16,7 +16,7 @@
 
 __all__ = ('mean_squared_error', 'mean_absolute_error',
            'mean_value_difference', 'standard_deviation_difference',
-           'range_difference', 'blurring', 'false_structures', 'ssim')
+           'range_difference', 'blurring', 'false_structures', 'ssim', 'psnr')
 
 
 def mean_squared_error(data, ground_truth, mask=None, normalized=False):
@@ -573,3 +573,62 @@ def smoothen(img):
         return 0.5 - ssim / 2
     else:
         return ssim
+
+
+def psnr(data, ground_truth, normalize=False):
+    """Return the Peak Signal-to-Noise Ratio.
+
+    Parameters
+    ----------
+    data : `FnBaseVector`
+        Input data or reconstruction.
+    ground_truth : `FnBaseVector`
+        Reference to compare ``data`` to.
+    normalize : bool
+        If true, normalizes ``data`` and ``ground_truth`` to have the same mean
+        and variance before comparison.
+
+    Returns
+    -------
+    psnr : float
+
+    Examples
+    --------
+    Compute the PSNR for two vectors:
+
+    >>> spc = odl.rn(5)
+    >>> data = spc.element([1, 1, 1, 1, 1])
+    >>> ground_truth = spc.element([1, 1, 1, 1, 2])
+    >>> result = psnr(data, ground_truth)
+    >>> print('{:.3f}'.format(result))
+    6.021
+
+    If data == ground_truth, the result is positive infinity:
+
+    >>> psnr(ground_truth, ground_truth)
+    inf
+
+    With ``normalize=True``, internal scaling and constant offsets are ignored:
+
+    >>> (psnr(data, ground_truth, normalize=True) ==
+    ...  psnr(data, 3 + 4 * ground_truth, normalize=True))
+    True
+    """
+    if normalize:
+        data = odl.util.zscore(data)
+        ground_truth = odl.util.zscore(ground_truth)
+
+    mse_result = mean_squared_error(data, ground_truth)
+    max_true = np.max(np.abs(ground_truth))
+
+    if mse_result == 0:
+        return np.inf
+    elif max_true == 0:
+        return -np.inf
+    else:
+        return 20 * np.log10(max_true) - 10 * np.log10(mse_result)
+
+if __name__ == '__main__':
+    # pylint: disable=wrong-import-position
+    from odl.util.testutils import run_doctests
+    run_doctests()

odl/contrib/fom/test/test_unsupervised.py

@@ -0,0 +1,57 @@
+# Copyright 2014-2017 The ODL contributors
+#
+# This file is part of ODL.
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License,
+# v. 2.0. If a copy of the MPL was not distributed with this file, You can
+# obtain one at https://mozilla.org/MPL/2.0/.
+
+"""Tests for unsupervised FoMs."""
+
+import numpy as np
+import pytest
+import odl
+import odl.contrib.fom
+
+
+def test_estimate_noise_std_constant_1d():
+    """Verify ``estimate_noise_std(0) == 0``"""
+    img = np.zeros(10)
+    result = odl.contrib.fom.estimate_noise_std(img)
+    assert pytest.approx(result) == 0.0
+
+
+def test_estimate_noise_std_normal_1d():
+    """Verify ``estimate_noise_std(N(0, 1)) == 1`` in 1d."""
+    img = np.random.randn(10)
+    result = odl.contrib.fom.estimate_noise_std(img)
+    expected = np.std(img)
+    assert pytest.approx(result, abs=0.2) == expected
+
+
+def test_estimate_noise_std_normal_2d():
+    """Verify ``estimate_noise_std(N(0, 1)) == 1`` in 2d."""
+    img = np.random.randn(10, 10)
+    result = odl.contrib.fom.estimate_noise_std(img)
+    expected = np.std(img)
+    assert pytest.approx(result, abs=0.2) == expected
+
+
+def test_estimate_noise_std_normal_4d():
+    """Verify ``estimate_noise_std(N(0, 1)) == 1`` in 4d."""
+    img = np.random.randn(5, 5, 5, 5)
+    result = odl.contrib.fom.estimate_noise_std(img)
+    expected = np.std(img)
+    assert pytest.approx(result, abs=0.2) == expected
+
+
+def test_estimate_noise_std_normal_large_1d():
+    """Verify ``estimate_noise_std(N(0, 1)) == 1`` with low error."""
+    img = np.random.randn(100000)
+    result = odl.contrib.fom.estimate_noise_std(img)
+    expected = np.std(img)
+    assert pytest.approx(result, abs=0.01) == expected
+
+
+if __name__ == '__main__':
+    pytest.main([str(__file__.replace('\\', '/')), '-v'])

odl/contrib/fom/unsupervised.py

@@ -0,0 +1,65 @@
+# Copyright 2014-2017 The ODL contributors
+#
+# This file is part of ODL.
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License,
+# v. 2.0. If a copy of the MPL was not distributed with this file, You can
+# obtain one at https://mozilla.org/MPL/2.0/.
+
+"""Figures of Merit (FOMs) for measuring image quality without a reference."""
+
+import numpy as np
+
+__all__ = ('estimate_noise_std',)
+
+
+def estimate_noise_std(img):
+    """Estimate standard deviation of noise in ``img``.
+
+    The algorithm, given in [Immerkaer1996], estimates the noise in an image
+    by
+
+    Parameters
+    ----------
+    img : array-like
+
+    Returns
+    -------
+    noise : float
+
+    Examples
+    --------
+    Create image with noise 1.0, verify result
+
+    >>> img = np.random.randn(10, 10)
+    >>> result = estimate_noise_std(img)  # should be about 1
+
+    Also works with higher dimensional arrays
+
+    >>> img = np.random.randn(3, 3, 3)
+    >>> result = estimate_noise_std(img)
+
+    References
+    ----------
+    [Immerkaer1996] Immerkaer, J. *Fast Noise Variance Estimation*.
+    Computer Vision and Image Understanding, 1996.
+    """
+    import scipy.signal
+    import functools
+    img = np.asarray(img, dtype='float')
+
+    M = functools.reduce(np.add.outer, [[-1, 2, -1]] * img.ndim)
+
+    convolved = scipy.signal.fftconvolve(img, M, mode='valid')
+    conv_var = np.sum(convolved ** 2)
+
+    scale = np.sum(np.square(M)) * convolved.size
+    sigma = np.sqrt(conv_var / scale)
+
+    return sigma
+
+
+if __name__ == '__main__':
+    # pylint: disable=wrong-import-position
+    from odl.util.testutils import run_doctests
+    run_doctests()

odl/contrib/tensorflow/README.md

@@ -0,0 +1,24 @@
+# Tensorflow
+
+This package contains ODL functionality related to [TensorFlow](https://www.tensorflow.org/).
+
+## Content
+
+* `TensorflowOperator` in [operator.py](operator.py) wraps a tensorflow expression into an ODL operator.
+  This allows using tensorflow neural networks as operators in ODL.
+* `as_tensorflow_layer` in [layer.py](layer.py) wraps an ODL operator into a tensorflow layer.
+  This allows using arbitrary ODL operators inside tensorflow neural networks.
+* `TensorflowSpace` in [space.py](space.py) is a `FnBase`-type space which uses tensorflow as a backend.
+
+## Example usage
+
+The [examples](examples) folder contains examples on how to use the above functionality.
+Specifically:
+
+* [tensorflow_layer_matrix.py](examples/tensorflow_layer_matrix.py) shows how an ODL `MatrixOperator` can be converted to a tensorflow layer.
+* [tensorflow_layer_productspace.py](examples/tensorflow_layer_productspace.py) shows how an ODL operator acting on `ProductSpace`s can be converted to a tensorflow layer.
+* [tensorflow_layer_ray_transform.py](examples/tensorflow_layer_ray_transform.py) shows how a `RayTransform` can be converted to a tensorflow layer.
+* [tensorflow_operator_matrix.py](examples/tensorflow_operator_matrix.py) shows how `tf.matmul` can be used as an ODL operator.
+* [tensorflow_tomography.py](examples/tensorflow_tomography.py) shows how tensorflow optimizers can be used with ODL operators to solve inverse problems.
+
+There are also some rudimentary tests in the [test](test) folder.

odl/contrib/tensorflow/__init__.py

@@ -0,0 +1,23 @@
+# Copyright 2014-2017 The ODL contributors
+#
+# This file is part of ODL.
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License,
+# v. 2.0. If a copy of the MPL was not distributed with this file, You can
+# obtain one at https://mozilla.org/MPL/2.0/.
+
+"""ODL interactions with tensorflow."""
+
+
+from __future__ import absolute_import
+
+__all__ = ()
+
+from .layer import *
+__all__ += layer.__all__
+
+from .space import *
+__all__ += space.__all__
+
+from .operator import *
+__all__ += operator.__all__

odl/contrib/tensorflow/examples/tensorflow_layer_matrix.py

@@ -0,0 +1,51 @@
+"""Example of how to convert an ODL operator to a tensorflow layer.
+
+In this example we take an ODL operator given by a `MatrixOperator` and
+convert it into a tensorflow layer that can be used inside any tensorflow
+computational graph.
+
+We also demonstrate that we can compute the "gradients" properly using the
+adjoint of the derivative.
+"""
+
+import tensorflow as tf
+import numpy as np
+import odl
+import odl.contrib.tensorflow
+
+sess = tf.InteractiveSession()
+
+# Define ODL operator
+matrix = np.array([[1, 2],
+                   [0, 0],
+                   [0, 1]], dtype='float32')
+odl_op = odl.MatrixOperator(matrix)
+
+# Define evaluation points
+x = [1., 2.]
+z = [1., 2., 3.]
+
+# Add empty axes for batch and channel
+x_tf = tf.constant(x)[None, ..., None]
+z_tf = tf.constant(z)[None, ..., None]
+
+# Create tensorflow layer from odl operator
+odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer(
+    odl_op, 'MatrixOperator')
+y_tf = odl_op_layer(x_tf)
+
+# Evaluate using tensorflow
+print('Tensorflow eval:                  ',
+      y_tf.eval().ravel())
+
+# Compare result with pure ODL
+print('ODL eval:                         ',
+      odl_op(x))
+
+# Evaluate the adjoint of the derivative, called gradient in tensorflow
+print('Tensorflow adjoint of derivative: ',
+      tf.gradients(y_tf, [x_tf], z_tf)[0].eval().ravel())
+
+# Compare result with pure ODL
+print('ODL adjoint of derivative:        ',
+      odl_op.derivative(x).adjoint(z))

odl/contrib/tensorflow/examples/tensorflow_layer_productspace.py

@@ -0,0 +1,46 @@
+"""Example of how to convert an ODL operator to a tensorflow layer.
+
+This example is similar to ``tensorflow_layer_matrix``, but demonstrates how
+to handle product-spaces.
+"""
+
+import tensorflow as tf
+import odl
+import odl.contrib.tensorflow
+
+sess = tf.InteractiveSession()
+
+# Define ODL operator
+space = odl.uniform_discr([0, 0], [1, 1], [10, 10], dtype='float32')
+
+odl_op = odl.Gradient(space)
+
+# Define evaluation points
+x = odl_op.domain.one()
+z = odl_op.range.one()
+
+# Add empty axes for batch and channel
+x_tf = tf.ones([1, 10, 10, 1])
+z_tf = tf.ones([1, 2, 10, 10, 1])
+
+# Create tensorflow layer from odl operator
+odl_op_layer = odl.contrib.tensorflow.as_tensorflow_layer(
+        odl_op, 'Gradient')
+y_tf = odl_op_layer(x_tf)
+
+# Evaluate using tensorflow
+print('Tensorflow eval:')
+print(y_tf.eval().squeeze())
+
+# Compare result with pure ODL
+print('ODL eval')
+print(odl_op(x))
+
+# Evaluate the adjoint of the derivative, called gradient in tensorflow
+scale = 1 / space.cell_volume
+print('Tensorflow adjoint of derivative (gradients):')
+print(tf.gradients(y_tf, [x_tf], z_tf)[0].eval().squeeze() * scale)
+
+# Compare result with pure ODL
+print('ODL adjoint of derivative:')
+print(odl_op.derivative(x).adjoint(z))