<commit_before># 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/>.

"""Tests for the Ray transform."""

# Imports for common Python 2/3 codebase
from __future__ import print_function, division, absolute_import
from future import standard_library
standard_library.install_aliases()

# External
import pytest
import numpy as np

# Internal
import odl
import odl.tomo as tomo
from odl.tomo.util.testutils import (skip_if_no_astra, skip_if_no_astra_cuda,
                                     skip_if_no_scikit)
from odl.util.testutils import almost_equal


# Find the valid projectors
projectors = [skip_if_no_astra('par2d astra_cpu uniform'),
              skip_if_no_astra('par2d astra_cpu nonuniform'),
              skip_if_no_astra('par2d astra_cpu random'),
              skip_if_no_astra('cone2d astra_cpu uniform'),
              skip_if_no_astra('cone2d astra_cpu nonuniform'),
              skip_if_no_astra('cone2d astra_cpu random'),
              skip_if_no_astra_cuda('par2d astra_cuda uniform'),
              skip_if_no_astra_cuda('par2d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('par2d astra_cuda random'),
              skip_if_no_astra_cuda('cone2d astra_cuda uniform'),
              skip_if_no_astra_cuda('cone2d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('cone2d astra_cuda random'),
              skip_if_no_astra_cuda('par3d astra_cuda uniform'),
              skip_if_no_astra_cuda('par3d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('par3d astra_cuda random'),
              skip_if_no_astra_cuda('cone3d astra_cuda uniform'),
              skip_if_no_astra_cuda('cone3d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('cone3d astra_cuda random'),
              skip_if_no_astra_cuda('helical astra_cuda uniform'),
              skip_if_no_scikit('par2d scikit uniform')]

projector_ids = ['geom={}, impl={}, angles={}'
                 ''.format(*p.args[1].split()) for p in projectors]


@pytest.fixture(scope="module", params=projectors, ids=projector_ids)
def projector(request):
    n_voxels = 100
    n_angles = 100
    n_pixels = 100

    geom, impl, angle = request.param.split()

    if angle == 'uniform':
        apart = odl.uniform_partition(0, 2 * np.pi, n_angles)
    elif angle == 'random':
        # Linearly spaced with random noise
        min_pt = 2 * (2.0 * np.pi) / n_angles
        max_pt = (2.0 * np.pi) - 2 * (2.0 * np.pi) / n_angles
        points = np.linspace(min_pt, max_pt, n_angles)
        points += np.random.rand(n_angles) * (max_pt - min_pt) / (5 * n_angles)
        agrid = odl.TensorGrid(points)
        apart = odl.RectPartition(odl.Interval(0, 2 * np.pi), agrid)
    elif angle == 'nonuniform':
        # Angles spaced quadratically
        min_pt = 2 * (2.0 * np.pi) / n_angles
        max_pt = (2.0 * np.pi) - 2 * (2.0 * np.pi) / n_angles
        points = np.linspace(min_pt ** 0.5, max_pt ** 0.5, n_angles) ** 2
        agrid = odl.TensorGrid(points)
        apart = odl.RectPartition(odl.Interval(0, 2 * np.pi), agrid)
    else:
        raise ValueError('angle not valid')

    if geom == 'par2d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 2, [20] * 2, [n_voxels] * 2,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition(-30, 30, n_pixels)
        geom = tomo.Parallel2dGeometry(apart, dpart)

        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'par3d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 3, [20] * 3, [n_voxels] * 3,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition([-30] * 2, [30] * 2, [n_pixels] * 2)
        geom = tomo.Parallel3dAxisGeometry(apart, dpart)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'cone2d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 2, [20] * 2, [n_voxels] * 2,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition(-30, 30, n_pixels)
        geom = tomo.FanFlatGeometry(apart, dpart, src_radius=200,
                                    det_radius=100)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'cone3d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 3, [20] * 3, [n_voxels] * 3,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition([-60] * 2, [60] * 2, [n_pixels] * 2)

        geom = tomo.CircularConeFlatGeometry(apart, dpart, src_radius=200,
                                             det_radius=100)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'helical':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20, -20, 0], [20, 20, 40],
                                       [n_voxels] * 3, dtype='float32')

        # overwrite angle
        apart = odl.uniform_partition(0, 8 * 2 * np.pi, n_angles)
        dpart = odl.uniform_partition([-30, -3], [30, 3], [n_pixels] * 2)
        geom = tomo.HelicalConeFlatGeometry(apart, dpart, pitch=5.0,
                                            src_radius=200, det_radius=100)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)
    else:
        raise ValueError('geom not valid')


def test_projector(projector):
    """Test discrete Ray transform forward projection."""

    # TODO: this needs to be improved
    # Accept 10% errors
    places = 1

    # Create Shepp-Logan phantom
    vol = projector.domain.one()

    # Calculate projection
    proj = projector(vol)

    # We expect maximum value to be along diagonal
    expected_max = projector.domain.partition.extent()[0] * np.sqrt(2)
    assert almost_equal(proj.ufunc.max(), expected_max, places=places)


def test_adjoint(projector):
    """Test discrete Ray transform backward projection."""

    # TODO: this needs to be improved<commit_msg>TST: add inplace test for tomo projectors, closes #266<commit_after># 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/>.

"""Tests for the Ray transform."""

# Imports for common Python 2/3 codebase
from __future__ import print_function, division, absolute_import
from future import standard_library
standard_library.install_aliases()

# External
import pytest
import numpy as np

# Internal
import odl
import odl.tomo as tomo
from odl.tomo.util.testutils import (skip_if_no_astra, skip_if_no_astra_cuda,
                                     skip_if_no_scikit)
from odl.util.testutils import almost_equal


# Find the valid projectors
projectors = [skip_if_no_astra('par2d astra_cpu uniform'),
              skip_if_no_astra('par2d astra_cpu nonuniform'),
              skip_if_no_astra('par2d astra_cpu random'),
              skip_if_no_astra('cone2d astra_cpu uniform'),
              skip_if_no_astra('cone2d astra_cpu nonuniform'),
              skip_if_no_astra('cone2d astra_cpu random'),
              skip_if_no_astra_cuda('par2d astra_cuda uniform'),
              skip_if_no_astra_cuda('par2d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('par2d astra_cuda random'),
              skip_if_no_astra_cuda('cone2d astra_cuda uniform'),
              skip_if_no_astra_cuda('cone2d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('cone2d astra_cuda random'),
              skip_if_no_astra_cuda('par3d astra_cuda uniform'),
              skip_if_no_astra_cuda('par3d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('par3d astra_cuda random'),
              skip_if_no_astra_cuda('cone3d astra_cuda uniform'),
              skip_if_no_astra_cuda('cone3d astra_cuda nonuniform'),
              skip_if_no_astra_cuda('cone3d astra_cuda random'),
              skip_if_no_astra_cuda('helical astra_cuda uniform'),
              skip_if_no_scikit('par2d scikit uniform')]

projector_ids = ['geom={}, impl={}, angles={}'
                 ''.format(*p.args[1].split()) for p in projectors]


@pytest.fixture(scope="module", params=projectors, ids=projector_ids)
def projector(request):
    n_voxels = 100
    n_angles = 100
    n_pixels = 100

    geom, impl, angle = request.param.split()

    if angle == 'uniform':
        apart = odl.uniform_partition(0, 2 * np.pi, n_angles)
    elif angle == 'random':
        # Linearly spaced with random noise
        min_pt = 2 * (2.0 * np.pi) / n_angles
        max_pt = (2.0 * np.pi) - 2 * (2.0 * np.pi) / n_angles
        points = np.linspace(min_pt, max_pt, n_angles)
        points += np.random.rand(n_angles) * (max_pt - min_pt) / (5 * n_angles)
        agrid = odl.TensorGrid(points)
        apart = odl.RectPartition(odl.Interval(0, 2 * np.pi), agrid)
    elif angle == 'nonuniform':
        # Angles spaced quadratically
        min_pt = 2 * (2.0 * np.pi) / n_angles
        max_pt = (2.0 * np.pi) - 2 * (2.0 * np.pi) / n_angles
        points = np.linspace(min_pt ** 0.5, max_pt ** 0.5, n_angles) ** 2
        agrid = odl.TensorGrid(points)
        apart = odl.RectPartition(odl.Interval(0, 2 * np.pi), agrid)
    else:
        raise ValueError('angle not valid')

    if geom == 'par2d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 2, [20] * 2, [n_voxels] * 2,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition(-30, 30, n_pixels)
        geom = tomo.Parallel2dGeometry(apart, dpart)

        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'par3d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 3, [20] * 3, [n_voxels] * 3,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition([-30] * 2, [30] * 2, [n_pixels] * 2)
        geom = tomo.Parallel3dAxisGeometry(apart, dpart)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'cone2d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 2, [20] * 2, [n_voxels] * 2,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition(-30, 30, n_pixels)
        geom = tomo.FanFlatGeometry(apart, dpart, src_radius=200,
                                    det_radius=100)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'cone3d':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20] * 3, [20] * 3, [n_voxels] * 3,
                                       dtype='float32')

        # Geometry
        dpart = odl.uniform_partition([-60] * 2, [60] * 2, [n_pixels] * 2)

        geom = tomo.CircularConeFlatGeometry(apart, dpart, src_radius=200,
                                             det_radius=100)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)

    elif geom == 'helical':
        # Discrete reconstruction space
        reco_space = odl.uniform_discr([-20, -20, 0], [20, 20, 40],
                                       [n_voxels] * 3, dtype='float32')

        # overwrite angle
        apart = odl.uniform_partition(0, 8 * 2 * np.pi, n_angles)
        dpart = odl.uniform_partition([-30, -3], [30, 3], [n_pixels] * 2)
        geom = tomo.HelicalConeFlatGeometry(apart, dpart, pitch=5.0,
                                            src_radius=200, det_radius=100)

        # Ray transform
        return tomo.RayTransform(reco_space, geom, impl=impl)
    else:
        raise ValueError('geom not valid')


@pytest.fixture(scope="module",
                params=[True, False],
                ids=[' inplace ', ' out of place '])
def inplace(request):
    return request.param


def test_projector(projector, inplace):
    """Test discrete Ray transform forward projection."""

    # TODO: this needs to be improved
    # Accept 10% errors
    places = 1

    # Create Shepp-Logan phantom
    vol = projector.domain.one()

    # Calculate projection
    if inplace:
        proj = projector.range.zero()
        projector(vol, out=proj)
    else:
        proj = projector(vol)

    # We expect maximum value to be along diagonal
    expected_max = projector.domain.partition.extent()[0] * np.sqrt(2)
    assert almost_equal(proj.ufunc.max(), expected_max, places=places)


def test_adjoint(projector):
    """Test discrete Ray transform backward projection."""

    # TODO: this needs to be improved
{'fraction_non_alphanumeric': 0.10437014329507378, 'fraction_numerical': 0.02602124474228274, 'mean_word_length': 2.434026927784578, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '114586', 'n_tokens_mistral': 4845, 'n_tokens_neox': 4514, 'n_words': 1129}