diff --git a/pyramid/numcore/__init__.py b/pyramid/numcore/__init__.py
deleted file mode 100644
index f4e72d8c44552843d3edf2a27b240c65c30329da..0000000000000000000000000000000000000000
--- a/pyramid/numcore/__init__.py
+++ /dev/null
@@ -1,15 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Package for numerical core routines which enable fast computations.
-
-Modules
--------
-phasemapper_core
- Provides numerical core routines for :class:`~.PhaseMapper` class.
-kernel_core
- Provides numerical core routines for the :class:`~.Kernel` class.
-
-Notes
------
-Packages are written in `pyx`-format for the use with :mod:`~.cython`.
-
-"""
diff --git a/pyramid/numcore/phasemapper_core.pyx b/pyramid/numcore/phasemapper_core.pyx
deleted file mode 100644
index 59be7c8ff4a84efec61e432707b54c81ae58c6f5..0000000000000000000000000000000000000000
--- a/pyramid/numcore/phasemapper_core.pyx
+++ /dev/null
@@ -1,156 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Numerical core routines for the :mod:`~.pyramid.phasemapper` module.
-
-Provides helper functions to speed up phase mapping calculations in the
-:mod:`~pyramid.phasemapper` module, by using C-speed for the for-loops and by omitting boundary
-and wraparound checks.
-
-"""
-
-
-import numpy as np
-import math
-
-cimport cython
-cimport numpy as np
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def phasemapper_real_convolve(
- unsigned int v_dim, unsigned int u_dim,
- float[:, :] v_phi, float[:, :] u_phi,
- float[:, :] v_mag, float[:, :] u_mag,
- float[:, :] phase, float threshold):
- '''Numerical core routine for the phase calculation using the real space approach.
-
- Parameters
- ----------
- v_dim, u_dim : int
- Dimensions of the projection along the two major axes.
- v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
- Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- v_mag, u_mag : :class:`~numpy.ndarray` (N=2)
- Magnetization components in `u`- and `v`-direction.
- phase : :class:`~numpy.ndarray` (N=2)
- Matrix in which the resulting magnetic phase map should be stored.
- threshold : float
- The `threshold` determines which pixels contribute to the magnetic phase.
-
- Returns
- -------
- None
-
- '''
- cdef unsigned int i, j, p, q, p_c, q_c
- cdef float u_m, v_m
- for j in range(v_dim):
- for i in range(u_dim):
- u_m = u_mag[j, i]
- v_m = v_mag[j, i]
- p_c = u_dim - 1 - i
- q_c = v_dim - 1 - j
- if abs(u_m) > threshold:
- for q in range(v_dim):
- for p in range(u_dim):
- phase[q, p] += u_m * u_phi[q_c + q, p_c + p]
- if abs(v_m) > threshold:
- for q in range(v_dim):
- for p in range(u_dim):
- phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def jac_dot_real_convolve(
- unsigned int v_dim, unsigned int u_dim,
- float[:, :] u_phi, float[:, :] v_phi,
- float[:] vector,
- float[:] result):
- '''Numerical core routine for the Jacobi matrix multiplication for the phase mapper.
-
- Parameters
- ----------
- v_dim, u_dim : int
- Dimensions of the projection along the two major axes.
- v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
- Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- vector : :class:`~numpy.ndarray` (N=1)
- Input vector which should be multiplied by the Jacobi-matrix.
- result : :class:`~numpy.ndarray` (N=1)
- Vector in which the result of the multiplication should be stored.
-
- Returns
- -------
- None
-
- Notes
- -----
- The strategy involves iterating over the magnetization first and over the kernel in an inner
- iteration loop.
-
- '''
- cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, r, u, v, size
- size = u_dim * v_dim # Number of pixels
- s = 0 # Current contributing pixel (numbered consecutively)
- # Iterate over all contributingh pixels:
- for j in range(v_dim):
- for i in range(u_dim):
- u_min = (u_dim-1) - i # u_max = u_min + u_dim
- v_min = (v_dim-1) - j # v_max = v_min + v_dim
- r = 0 # Current result component / affected pixel (numbered consecutively)
- # Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
- for v in range(v_min, v_min + v_dim):
- for u in range(u_min, u_min + u_dim):
- result[r] += vector[s] * u_phi[v, u]
- result[r] -= vector[s+size] * v_phi[v, u]
- r += 1
- s += 1
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def jac_T_dot_real_convolve(
- unsigned int v_dim, unsigned int u_dim,
- float[:, :] u_phi, float[:, :] v_phi,
- float[:] vector,
- float[:] result):
- '''Core routine for the transposed Jacobi multiplication for the phase mapper.
-
- Parameters
- ----------
- v_dim, u_dim : int
- Dimensions of the projection along the two major axes.
- v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
- Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- vector : :class:`~numpy.ndarray` (N=1)
- Input vector which should be multiplied by the transposed Jacobi-matrix.
- result : :class:`~numpy.ndarray` (N=1)
- Vector in which the result of the multiplication should be stored.
-
- Returns
- -------
- None
-
- Notes
- -----
- The strategy involves iterating over the magnetization first and over the kernel in an inner
- iteration loop.
-
- '''
- cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, r, u, v, size
- size = u_dim * v_dim # Number of pixels
- r = 0 # Current result component / contributing pixel (numbered consecutively)
- # Iterate over all contributingh pixels:
- for j in range(v_dim):
- for i in range(u_dim):
- u_min = (u_dim-1) - i # u_max = u_min + u_dim
- v_min = (v_dim-1) - j # v_max = v_min + v_dim
- s = 0 # Current affected pixel (numbered consecutively)
- # Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
- for v in range(v_min, v_min + v_dim):
- for u in range(u_min, u_min + u_dim):
- result[r] += vector[s] * u_phi[v, u]
- result[r+size] -= vector[s] * v_phi[v, u]
- s += 1
- r += 1
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index c88e720268364ed204d5d3cec37739c68604992c..bded43f7a7334e0479aa5902fa6e608239cb8b20 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -11,7 +11,6 @@ import abc
import logging
import numpy as np
-import pyramid.numcore.phasemapper_core as nc
from numpy import pi
from pyramid import fft
@@ -20,7 +19,7 @@ from pyramid.fielddata import VectorData, ScalarData
from pyramid.phasemap import PhaseMap
from pyramid.projector import RotTiltProjector, XTiltProjector, YTiltProjector, SimpleProjector
-__all__ = ['PhaseMapperRDFC', 'PhaseMapperRDRC', 'PhaseMapperFDFC', 'PhaseMapperMIP', 'pm']
+__all__ = ['PhaseMapperRDFC', 'PhaseMapperFDFC', 'PhaseMapperMIP', 'pm']
_log = logging.getLogger(__name__)
PHI_0 = 2067.83 # magnetic flux in T*nm²
@@ -193,159 +192,6 @@ class PhaseMapperRDFC(PhaseMapper):
return result
-class PhaseMapperRDRC(PhaseMapper):
- """Class representing a phase mapping strategy using real space discretization.
-
- The :class:`~.PMReal` class represents a phase mapping strategy involving discretization in
- real space. It directly takes :class:`~.VectorData` objects and returns :class:`~.PhaseMap`
- objects.
-
- Attributes
- ----------
- kernel : :class:`~pyramid.Kernel`
- Convolution kernel, representing the phase contribution of one single magnetized pixel.
- threshold : float, optional
- Threshold determining for which voxels the phase contribution will be calculated. The
- default is 0, which means that all voxels with non-zero magnetization will contribute.
- Should be above noise level.
- numcore : boolean, optional
- Boolean choosing if Cython enhanced routines from the :mod:`~.pyramid.numcore` module
- should be used. Default is True.
- m: int
- Size of the image space.
- n: int
- Size of the input space.
-
- """
-
- _log = logging.getLogger(__name__ + '.PhaseMapperRDRC')
-
- def __init__(self, kernel, threshold=0, numcore=True):
- self._log.debug('Calling __init__')
- self.kernel = kernel
- self.threshold = threshold
- self.numcore = numcore
- self.m = np.prod(kernel.dim_uv)
- self.n = 2 * self.m
- self._log.debug('Created ' + str(self))
-
- def __repr__(self):
- self._log.debug('Calling __repr__')
- return '%s(kernel=%r, threshold=%r, numcore=%r)' % \
- (self.__class__, self.kernel, self.threshold, self.numcore)
-
- def __str__(self):
- self._log.debug('Calling __str__')
- return 'PhaseMapperRDRC(kernel=%s, threshold=%s, numcore=%s)' % \
- (self.kernel, self.threshold, self.numcore)
-
- def __call__(self, mag_data):
- self._log.debug('Calling __call__')
- dim_uv = self.kernel.dim_uv
- assert isinstance(mag_data, VectorData), 'Only VectorData objects can be mapped!'
- assert mag_data.a == self.kernel.a, 'Grid spacing has to match!'
- assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
- assert mag_data.dim[1:3] == dim_uv, 'Dimensions do not match!'
- # Process input parameters:
- u_mag, v_mag = mag_data.field[0:2, 0, ...]
- # Get kernel (lookup-tables for the phase of one pixel):
- u_phi = self.kernel.u
- v_phi = self.kernel.v
- # Calculation of the phase:
- phase = np.zeros(dim_uv, dtype=np.float32)
- if self.numcore:
- nc.phasemapper_real_convolve(dim_uv[0], dim_uv[1], u_phi, v_phi,
- v_mag, u_mag, phase, self.threshold)
- else:
- for j in range(dim_uv[0]):
- for i in range(dim_uv[1]):
- v_phase = v_phi[dim_uv[0] - 1 - j:(2 * dim_uv[0] - 1) - j,
- dim_uv[1] - 1 - i:(2 * dim_uv[1] - 1) - i]
- if abs(v_mag[j, i]) > self.threshold:
- phase += u_mag[j, i] * v_phase
- u_phase = u_phi[dim_uv[0] - 1 - j:(2 * dim_uv[0] - 1) - j,
- dim_uv[1] - 1 - i:(2 * dim_uv[1] - 1) - i]
- if abs(u_mag[j, i]) > self.threshold:
- phase -= v_mag[j, i] * u_phase
- # Return the phase:
- return PhaseMap(mag_data.a, phase)
-
- def jac_dot(self, vector):
- """Calculate the product of the Jacobi matrix with a given `vector`.
-
- Parameters
- ----------
- vector : :class:`~numpy.ndarray` (N=1)
- Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
- (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
- ``N**2`` elements to the `v`-component of the magnetization.
-
- Returns
- -------
- result : :class:`~numpy.ndarray` (N=1)
- Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
-
- """
- assert len(vector) == self.n, \
- 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
- v_dim, u_dim = self.kernel.dim_uv
- result = np.zeros(self.m, dtype=np.float32)
- if self.numcore:
- if vector.dtype != np.float32:
- vector = vector.astype(np.float32)
- nc.jac_dot_real_convolve(v_dim, u_dim, self.kernel.v, self.kernel.u, vector, result)
- else:
- # Iterate over all contributing pixels (numbered consecutively)
- for s in range(self.m): # column-wise (two columns at a time, u- and v-component)
- i = s % u_dim # u-coordinate of current contributing pixel
- j = int(s / u_dim) # v-coordinate of current ccontributing pixel
- u_min = (u_dim - 1) - i # u_dim-1: center of the kernel
- u_max = (2 * u_dim - 1) - i # = u_min + u_dim
- v_min = (v_dim - 1) - j # v_dim-1: center of the kernel
- v_max = (2 * v_dim - 1) - j # = v_min + v_dim
- result += vector[s] * self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1)
- result -= vector[s + self.m] * self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1)
- return result
-
- def jac_T_dot(self, vector):
- """Calculate the product of the transposed Jacobi matrix with a given `vector`.
-
- Parameters
- ----------
- vector : :class:`~numpy.ndarray` (N=1)
- Vector with ``N**2`` entries which represents a matrix with dimensions like a scalar
- phasemap.
-
- Returns
- -------
- result : :class:`~numpy.ndarray` (N=1)
- Product of the transposed Jacobi matrix (which is not explicitely calculated) with
- the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
-
- """
- assert len(vector) == self.m, \
- 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
- v_dim, u_dim = self.kernel.dim_uv
- result = np.zeros(self.n, dtype=np.float32)
- if self.numcore:
- if vector.dtype != np.float32:
- vector = vector.astype(np.float32)
- nc.jac_T_dot_real_convolve(v_dim, u_dim, self.kernel.v, self.kernel.u, vector, result)
- else:
- # Iterate over all contributing pixels (numbered consecutively):
- for s in range(self.m): # row-wise (two rows at a time, u- and v-component)
- i = s % u_dim # u-coordinate of current contributing pixel
- j = int(s / u_dim) # v-coordinate of current contributing pixel
- u_min = (u_dim - 1) - i # u_dim-1: center of the kernel
- u_max = (2 * u_dim - 1) - i # = u_min + u_dim
- v_min = (v_dim - 1) - j # v_dim-1: center of the kernel
- v_max = (2 * v_dim - 1) - j # = v_min + v_dim
- result[s] = np.sum(vector * self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1))
- result[s + self.m] = np.sum(vector *
- -self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1))
- return result
-
-
class PhaseMapperFDFC(PhaseMapper):
"""Class representing a phase mapping strategy using a discretization in Fourier space.
diff --git a/pyramid/tests/test_phasemapper.py b/pyramid/tests/test_phasemapper.py
index 9666546ae2f91dae658da13967ad06fcebade29a..d16745877ea2e785b05333d281c4d68b38775844 100644
--- a/pyramid/tests/test_phasemapper.py
+++ b/pyramid/tests/test_phasemapper.py
@@ -10,8 +10,7 @@ from numpy.testing import assert_allclose
from pyramid.kernel import Kernel
from pyramid.fielddata import VectorData, ScalarData
from pyramid.phasemap import PhaseMap
-from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperRDRC, PhaseMapperFDFC
-from pyramid.phasemapper import PhaseMapperMIP, pm
+from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperMIP, pm
class TestCasePhaseMapperRDFC(unittest.TestCase):
@@ -52,44 +51,6 @@ class TestCasePhaseMapperRDFC(unittest.TestCase):
err_msg='Unexpected behaviour in the the transposed jacobi matrix!')
-class TestCasePhaseMapperRDRC(unittest.TestCase):
- def setUp(self):
- self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
- self.mag_proj = VectorData.load_from_hdf5(os.path.join(self.path, 'mag_proj.hdf5'))
- self.mapper = PhaseMapperRDRC(Kernel(self.mag_proj.a, self.mag_proj.dim[1:]))
-
- def tearDown(self):
- self.path = None
- self.mag_proj = None
- self.mapper = None
-
- def test_PhaseMapperRDRC_call(self):
- phase_ref = PhaseMap.load_from_hdf5(os.path.join(self.path, 'phase_map.hdf5'))
- phase_map = self.mapper(self.mag_proj)
- assert_allclose(phase_map.phase, phase_ref.phase, atol=1E-7,
- err_msg='Unexpected behavior in __call__()!')
- assert_allclose(phase_map.a, phase_ref.a, err_msg='Unexpected behavior in __call__()!')
-
- def test_PhaseMapperRDRC_jac_dot(self):
- phase = self.mapper(self.mag_proj).phase
- mag_proj_vec = self.mag_proj.field[:2, ...].flatten()
- phase_jac = self.mapper.jac_dot(mag_proj_vec).reshape(self.mapper.kernel.dim_uv)
- assert_allclose(phase, phase_jac, atol=1E-7,
- err_msg='Inconsistency between __call__() and jac_dot()!')
- n = self.mapper.n
- jac = np.array([self.mapper.jac_dot(np.eye(n)[:, i]) for i in range(n)]).T
- jac_ref = np.load(os.path.join(self.path, 'jac.npy'))
- assert_allclose(jac, jac_ref, atol=1E-7,
- err_msg='Unexpected behaviour in the the jacobi matrix!')
-
- def test_PhaseMapperRDRC_jac_T_dot(self):
- m = self.mapper.m
- jac_T = np.array([self.mapper.jac_T_dot(np.eye(m)[:, i]) for i in range(m)]).T
- jac_T_ref = np.load(os.path.join(self.path, 'jac.npy')).T
- assert_allclose(jac_T, jac_T_ref, atol=1E-7,
- err_msg='Unexpected behaviour in the the transposed jacobi matrix!')
-
-
class TestCasePhaseMapperFDFCpad0(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
diff --git a/setup.py b/setup.py
index 73510a518183fa72a87bd44b6c48b8198333c2b5..da4440d35c1706967cccb4d58b24387a68b29a53 100644
--- a/setup.py
+++ b/setup.py
@@ -147,11 +147,8 @@ setup(name=DISTNAME,
packages=find_packages(exclude=['tests']),
include_dirs=[numpy.get_include()],
requires=['numpy', 'scipy', 'matplotlib', 'Pillow',
- 'mayavi', 'pyfftw', 'hyperspy', 'Cython', 'nose'],
+ 'mayavi', 'pyfftw', 'hyperspy', 'nose'],
scripts=get_files('scripts'),
test_suite='nose.collector',
- cmdclass={'build_ext': build_ext, 'build': build},
- ext_package='pyramid/numcore',
- ext_modules=[Extension('phasemapper_core', ['pyramid/numcore/phasemapper_core.pyx'],
- include_dirs=[numpy.get_include()])])
+ cmdclass={'build_ext': build_ext, 'build': build})
print('-------------------------------------------------------------------------------\n')