# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module executes several forward models to calculate the magnetic or electric phase map from
a given projection of a 3-dimensional magnetic distribution (see :mod:`~pyramid.projector`).
For the magnetic phase map, an approach using real space and one using Fourier space is provided.
The electrostatic contribution is calculated by using the assumption of a mean inner potential."""
import abc
import logging
import numpy as np
import pyramid.numcore.phasemapper_core as nc
from numpy import pi
from pyramid import fft
from pyramid.kernel import Kernel
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import RotTiltProjector, XTiltProjector, YTiltProjector, SimpleProjector
__all__ = ['PhaseMapperRDFC', 'PhaseMapperRDRC', 'PhaseMapperFDFC', 'pm']
_log = logging.getLogger(__name__)
PHI_0 = 2067.83 # magnetic flux in T*nm²
H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C
C = 2.998E8 # speed of light in m/s
class PhaseMapper(object):
"""Abstract base class for the phase calculation from a 2-dimensional distribution.
The :class:`~.PhaseMapper-` class represents a strategy for the phasemapping of a
2-dimensional magnetic distribution with two components onto a scalar phase map.
:class:`~.Kernel` is an abstract base class and provides a unified interface which should be
subclassed with custom :func:`__init__` and :func:`__call__` functions. Concrete subclasses
can be called as a function and take a :class:`~.MagData` object as input and return a
:class:`~.PhaseMap` object.
"""
__metaclass__ = abc.ABCMeta
_log = logging.getLogger(__name__ + '.PhaseMapper')
@abc.abstractmethod
def __call__(self, mag_data):
raise NotImplementedError()
@abc.abstractmethod
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.
"""
raise NotImplementedError()
@abc.abstractmethod
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.
"""
raise NotImplementedError()
class PhaseMapperRDFC(PhaseMapper):
"""Class representing a phase mapping strategy using real space discretization and Fourier
space convolution.
The :class:`~.PMConvolve` class represents a phase mapping strategy involving discretization in
real space. It utilizes the convolution in Fourier space, directly takes :class:`~.MagData`
objects and returns :class:`~.PhaseMap` objects.
Attributes
----------
kernel : :class:`~pyramid.Kernel`
Convolution kernel, representing the phase contribution of one single magnetized pixel.
m: int
Size of the image space.
n: int
Size of the input space.
"""
_log = logging.getLogger(__name__ + '.PhaseMapperRDFC')
def __init__(self, kernel):
self._log.debug('Calling __init__')
self.kernel = kernel
self.m = np.prod(kernel.dim_uv)
self.n = 2 * self.m
self.u_mag = fft.zeros(kernel.dim_pad, dtype=fft.FLOAT)
self.v_mag = fft.zeros(kernel.dim_pad, dtype=fft.FLOAT)
self.mag_adj = fft.zeros(kernel.dim_pad, dtype=fft.FLOAT)
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(kernel=%r)' % (self.__class__, self.kernel)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMapperRDFC(kernel=%s)' % self.kernel
def __call__(self, mag_data):
assert isinstance(mag_data, MagData), 'Only MagData 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] == self.kernel.dim_uv, 'Dimensions do not match!'
# Process input parameters:
self.u_mag[self.kernel.slice_mag] = mag_data.magnitude[0, 0, ...] # u-component
self.v_mag[self.kernel.slice_mag] = mag_data.magnitude[1, 0, ...] # v-component
return PhaseMap(mag_data.a, self._convolve())
def _convolve(self):
# Fourier transform the projected magnetisation:
self.u_mag_fft = fft.rfftn(self.u_mag)
self.v_mag_fft = fft.rfftn(self.v_mag)
# Convolve the magnetization with the kernel in Fourier space:
self.phase_fft = self.u_mag_fft * self.kernel.v_fft - self.v_mag_fft * self.kernel.u_fft
# Return the result:
return fft.irfftn(self.phase_fft)[self.kernel.slice_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)
self.u_mag[self.kernel.slice_mag], self.v_mag[self.kernel.slice_mag] = \
np.reshape(vector, (2,) + self.kernel.dim_uv)
return np.ravel(self._convolve())
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)
self.mag_adj[self.kernel.slice_mag] = vector.reshape(self.kernel.dim_uv)
mag_adj_fft = fft.irfftn_adj(self.mag_adj)
u_phase_adj_fft = mag_adj_fft * self.kernel.v_fft
v_phase_adj_fft = mag_adj_fft * -self.kernel.u_fft
u_phase_adj = fft.rfftn_adj(u_phase_adj_fft)[self.kernel.slice_phase]
v_phase_adj = fft.rfftn_adj(v_phase_adj_fft)[self.kernel.slice_phase]
result = np.concatenate((-u_phase_adj.flatten(), -v_phase_adj.flatten()))
# TODO: Why minus?
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:`~.MagData` 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, MagData), 'Only MagData 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.magnitude[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.
The :class:`~.PMFourier` class represents a phase mapping strategy involving discretization in
Fourier space. It utilizes the Fourier transforms, which are inherently calculated in the
:class:`~.Kernel` class and directly takes :class:`~.MagData` objects and returns
:class:`~.PhaseMap` objects.
Attributes
----------
a : float
The grid spacing in nm.
dim_uv : tuple of int (N=2)
Dimensions of the 2-dimensional projected magnetization grid for the kernel setup.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
padding : int, optional
Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
m: int
Size of the image space.
n: int
Size of the input space.
"""
_log = logging.getLogger(__name__ + '.PhaseMapperFDFC')
def __init__(self, a, dim_uv, b_0=1, padding=0):
self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.b_0 = b_0
self.padding = padding
self.m = np.prod(dim_uv)
self.n = 2 * self.m
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, b_0=%r, padding=%r)' % \
(self.__class__, self.a, self.dim_uv, self.b_0, self.padding)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMapperFDFC(a=%s, dim_uv=%s, b_0=%s, padding=%s)' % \
(self.a, self.dim_uv, self.b_0, self.padding)
def __call__(self, mag_data):
self._log.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
v_dim, u_dim = self.dim_uv
u_mag, v_mag = mag_data.magnitude[0:2, 0, ...]
# Create zero padded matrices:
u_pad = int(u_dim / 2 * self.padding)
v_pad = int(v_dim / 2 * self.padding)
u_mag_pad = np.pad(u_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant')
v_mag_pad = np.pad(v_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant')
# Fourier transform of the two components:
u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_pad), axes=0)
v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_pad), axes=0)
# Calculate the Fourier transform of the phase:
f_nyq = 0.5 / self.a # nyquist frequency
f_u = np.linspace(0, f_nyq, u_mag_fft.shape[1])
f_v = np.linspace(-f_nyq, f_nyq, u_mag_fft.shape[0], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
coeff = - (1j * self.b_0 * self.a) / (2 * PHI_0) # Minus because of negative z-direction
phase_fft = coeff * (u_mag_fft * f_vv - v_mag_fft * f_uu) / (f_uu ** 2 + f_vv ** 2 + 1e-30)
# Transform to real space and revert padding:
phase_pad = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
phase = phase_pad[v_pad:v_pad + v_dim, u_pad:u_pad + u_dim]
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.
"""
self._log.debug('Calling jac_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
mag_proj = MagData(self.a, np.zeros((3, 1) + self.dim_uv, dtype=np.float32))
magnitude_proj = np.reshape(vector, (2,) + self.dim_uv)
mag_proj.magnitude[:2, 0, ...] = magnitude_proj
return self(mag_proj).phase_vec
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.
"""
raise NotImplementedError()
class PhaseMapperMIP(PhaseMapper):
"""Class representing a phase mapping strategy for the electrostatic contribution.
The :class:`~.PhaseMapperMIP` class represents a phase mapping strategy for the electrostatic
contribution to the electron phase shift which results e.g. from the mean inner potential in
certain samples and which is sensitive to properties of the electron microscope. It directly
takes :class:`~.MagData` objects and returns :class:`~.PhaseMap` objects.
Attributes
----------
a : float
The grid spacing in nm.
dim_uv : tuple of int (N=2)
Dimensions of the 2-dimensional projected magnetization grid for the kernel setup.
v_0 : float, optional
The mean inner potential of the specimen in V. The default is 1.
v_acc : float, optional
The acceleration voltage of the electron microscope in V. The default is 30000.
threshold : float, optional
Threshold for the recognition of the specimen location. The default is 0, which means that
all voxels with non-zero magnetization will contribute. Should be above noise level.
m: int
Size of the image space.
n: int
Size of the input space.
"""
_log = logging.getLogger(__name__ + '.PhaseMapperMIP')
def __init__(self, a, dim_uv, v_0=1, v_acc=30000, threshold=0):
self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.v_0 = v_0
self.v_acc = v_acc
self.threshold = threshold
self.m = np.prod(self.dim_uv)
self.n = np.prod(self.dim_uv)
# Coefficient calculation:
lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E * v_acc / (2 * M_E * C ** 2)))
C_e = 2 * pi * Q_E / lam * (Q_E * v_acc + M_E * C ** 2) / (
Q_E * v_acc * (Q_E * v_acc + 2 * M_E * C ** 2))
self.coeff = v_0 * C_e * a * 1E-9
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, v_0=%r, v_acc=%r, threshold=%r)' % \
(self.__class__, self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMapperMIP(a=%s, dim_uv=%s, v_0=%s, v_acc=%s, threshold=%s)' % \
(self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
def __call__(self, mag_data):
self._log.debug('Calling __call__')
raise NotImplementedError() # TODO: Implement right!
# assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
# assert mag_data.a == self.a, 'Grid spacing has to match!'
# assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
# assert mag_data.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
# return self.coeff * mag_data.get_mask(self.threshold)[0, ...].reshape(self.dim_uv)
# # Calculate mask:
# mask = mag_data.get_mask(self.threshold)
# # Project and calculate phase:
# # TODO: check if projector manages scalar multiplication (problem with MagData) and fix!
# projection = self.projector(mask.reshape(-1)).reshape(self.projector.dim_uv)
# phase = self.coeff * projection
# 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.
"""
raise NotImplementedError()
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.
"""
raise NotImplementedError()
def pm(mag_data, mode='z', b_0=1, **kwargs):
"""Convenience function for fast phase mapping.
Parameters
----------
mag_data : :class:`~.MagData`
A :class:`~.MagData` object, from which the projected phase map should be calculated.
mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
is used for the projection. Default is a simple projection along the `z`-direction.
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
**kwargs : additional arguments
Additional arguments like `dim_uv`, 'tilt' or 'rotation', which are passed to the
projector-constructor, defined by the `mode`.
Returns
-------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
The calculated phase map as a :class:`~.PhaseMap` object.
"""
_log.debug('Calling pm')
# Determine projection mode:
if mode == 'rot-tilt':
projector = RotTiltProjector(mag_data.dim, **kwargs)
elif mode == 'x-tilt':
projector = XTiltProjector(mag_data.dim, **kwargs)
elif mode == 'y-tilt':
projector = YTiltProjector(mag_data.dim, **kwargs)
elif mode in ['x', 'y', 'z']:
projector = SimpleProjector(mag_data.dim, axis=mode, **kwargs)
else:
raise ValueError("Invalid mode (use 'x', 'y', 'z', 'x-tilt', 'y-tilt' or 'rot-tilt')")
# Project:
mag_proj = projector(mag_data)
# Set up phasemapper and map phase:
phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv, b_0=b_0))
phase_map = phasemapper(mag_proj)
# Get mask from magdata:
phase_map.mask = mag_proj.get_mask()[0, ...]
# Return phase:
return phase_map