Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • empyre/empyre
  • weber/empyre
  • wessels/empyre
  • bryan/empyre
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 6193 deletions
pyramid/analytic.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Create phase maps for magnetic distributions with analytic solutions.
This module provides methods for the calculation of the magnetic phase for simple geometries for
which the analytic solutions are known. These can be used for comparison with the phase
calculated by the functions from the :mod:`~pyramid.phasemapper` module.
"""
import logging
import numpy as np
from numpy import pi
from pyramid.phasemap import PhaseMap
__all__ = ['phase_mag_slab', 'phase_mag_slab', 'phase_mag_sphere', 'phase_mag_vortex']
_log = logging.getLogger(__name__)
PHI_0 = 2067.83 # magnetic flux in T*nm²
def phase_mag_slab(dim, a, phi, center, width, b_0=1):
"""Calculate the analytic magnetic phase for a homogeneously magnetized slab.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the slab in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the slab in pixel coordinates `(z, y, x)`.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phasemap : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
"""
_log.debug('Calling phase_mag_slab')
# Function for the phase:
def _phi_mag(x, y):
def _F_0(x, y):
A = np.log(x ** 2 + y ** 2 + 1E-30)
B = np.arctan(x / (y + 1E-30))
return x * A - 2 * x + 2 * y * B
return coeff * Lz * (- np.cos(phi) * (_F_0(x - x0 - Lx / 2, y - y0 - Ly / 2) -
_F_0(x - x0 + Lx / 2, y - y0 - Ly / 2) -
_F_0(x - x0 - Lx / 2, y - y0 + Ly / 2) +
_F_0(x - x0 + Lx / 2, y - y0 + Ly / 2))
+ np.sin(phi) * (_F_0(y - y0 - Ly / 2, x - x0 - Lx / 2) -
_F_0(y - y0 + Ly / 2, x - x0 - Lx / 2) -
_F_0(y - y0 - Ly / 2, x - x0 + Lx / 2) +
_F_0(y - y0 + Ly / 2, x - x0 + Lx / 2)))
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * center[1] # y0, x0 define the center of a pixel,
x0 = a * center[2] # hence: (cellindex + 0.5) * grid spacing
Lz, Ly, Lx = a * width[0], a * width[1], a * width[2]
coeff = - b_0 / (4 * PHI_0) # Minus because of negative z-direction
# Create grid:
x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, _phi_mag(xx, yy))
def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
"""Calculate the analytic magnetic phase for a homogeneously magnetized disc.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phasemap : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
"""
_log.debug('Calling phase_mag_disc')
# Function for the phase:
def _phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
result = coeff * Lz * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r <= R, 1, (R / (r + 1E-30)) ** 2)
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * center[1]
x0 = a * center[2]
Lz = a * height
R = a * radius
coeff = pi * b_0 / (2 * PHI_0) # Minus is gone because of negative z-direction
# Create grid:
x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, _phi_mag(xx, yy))
def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
"""Calculate the analytic magnetic phase for a homogeneously magnetized sphere.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the sphere in pixel coordinates `(z, y, x)`.
radius : float
The radius of the sphere in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phasemap : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
"""
_log.debug('Calling phase_mag_sphere')
# Function for the phase:
def _phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
result = coeff * R ** 3 / (r + 1E-30) ** 2 * (
(y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * center[1]
x0 = a * center[2]
R = a * radius
coeff = 2. / 3. * pi * b_0 / PHI_0 # Minus is gone because of negative z-direction
# Create grid:
x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, _phi_mag(xx, yy))
def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
"""Calculate the analytic magnetic phase for a vortex state disc.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`, which is also the vortex center.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phasemap : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
"""
_log.debug('Calling phase_mag_vortex')
# Function for the phase:
def _phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
result = coeff * np.where(r <= R, r - R, 0)
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * center[1]
x0 = a * center[2]
Lz = a * height
R = a * radius
coeff = - pi * b_0 * Lz / PHI_0 # Minus because of negative z-direction
# Create grid:
x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, _phi_mag(xx, yy))
pyramid/colors.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides a custom :class:`~.HLSTriadicColormap` colormap class which has a few
additional functions and can encode three-dimensional directions."""
import logging
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from matplotlib import colors
import colorsys
__all__ = ['TransparentColormap', 'HLSTriadicColormap', 'HLSTetradicColormap', 'hls_triadic_cmap',
'hls_tetradic_cmap', 'plot_colorwheel', 'rgb_from_hls',
'hls_from_vector', 'rgb_from_vector']
_log = logging.getLogger(__name__)
class TransparentColormap(colors.LinearSegmentedColormap):
"""Colormap subclass for including transparency.
This class is a subclass of the :class:`~matplotlib.pyplot.colors.LinearSegmentedColormap`
class with integrated support for transparency. The colormap is unicolor and varies only in
transparency.
Attributes
----------
r: float, optional
Intensity of red in the colormap. Has to be between 0. and 1.
g: float, optional
Intensity of green in the colormap. Has to be between 0. and 1.
b: float, optional
Intensity of blue in the colormap. Has to be between 0. and 1.
alpha_range : list (N=2) of float, optional
Start and end alpha value. Has to be between 0. and 1.
"""
_log = logging.getLogger(__name__ + '.TransparentColormap')
def __init__(self, r=1., g=0., b=0., alpha_range=None):
self._log.debug('Calling __init__')
if alpha_range is None:
alpha_range = [0., 1.]
red = [(0., 0., r), (1., r, 1.)]
green = [(0., 0., g), (1., g, 1.)]
blue = [(0., 0., b), (1., b, 1.)]
alpha = [(0., 0., alpha_range[0]), (1., alpha_range[1], 1.)]
cdict = {'red': red, 'green': green, 'blue': blue, 'alpha': alpha}
super().__init__('transparent', cdict, N=256)
self._log.debug('Created ' + str(self))
class HLSTriadicColormap(colors.ListedColormap):
"""Colormap subclass for encoding directions with colors.
This class is a subclass of the :class:`~matplotlib.pyplot.colors.ListedColormap`
class. The class follows the HSL ('hue', 'saturation', 'lightness') 'Double Hexcone' Model
with the saturation always set to 1 (moving on the surface of the color
cylinder) with a luminance of 0.5 (full color). The three prime colors (`rgb`) are spaced
equidistant with 120° space in between, according to a triadic arrangement.
"""
_log = logging.getLogger(__name__ + '.HLSTriadicColormap')
def __init__(self):
self._log.debug('Calling __init__')
h = np.linspace(0, 1, 256)
l = 0.5 * np.ones_like(h)
s = np.ones_like(h)
rgb = rgb_from_hls(h, l, s)
colors = [(ri / 255, gi / 255, bi / 255) for ri, gi, bi in rgb]
super().__init__(colors, 'hlscm', N=256)
self._log.debug('Created ' + str(self))
class HLSTetradicColormap(colors.LinearSegmentedColormap):
"""Colormap subclass for encoding directions with colors.
This class is a subclass of the :class:`~matplotlib.pyplot.colors.LinearSegmentedColormap`
class. The class follows the HSL ('hue', 'saturation', 'lightness') 'Double
Hexcone' Model with the saturation always set to 1 (moving on the surface of the color
cylinder) with a luminance of 0.5 (full color). The colors follow a tetradic arrangement with
four colors (red, green, blue and yellow) arranged with 90° spacing in between.
"""
_log = logging.getLogger(__name__ + '.HLSTetradicColormap')
CDICT = {'red': [(0.00, 1.0, 1.0),
(0.25, 0.0, 0.0),
(0.50, 0.0, 0.0),
(0.75, 1.0, 1.0),
(1.00, 1.0, 1.0)],
'green': [(0.00, 0.0, 0.0),
(0.25, 0.0, 0.0),
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 0.0)],
'blue': [(0.00, 0.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 0.0, 0.0)]}
def __init__(self):
self._log.debug('Calling __init__')
super().__init__('holo', self.CDICT, N=256)
self._log.debug('Created ' + str(self))
def plot_colorwheel(mode='triadic'):
"""Display a color wheel to illustrate the color coding of vector gradient directions.
Parameters
----------
mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme of the color wheel. The default is the
standard HLS scheme with the three primary colors (red, blue, green) spaced with 120°
between them (triadic arrangement). The other option is a tetradic scheme commonly
used in holography where red, green, blue and yellow are used with 90° spacing. In both
cases, the saturation decreases to the center of the circle.
Returns
-------
None
"""
_log.debug('Calling display_color_wheel')
# Construct the colorwheel:
yy, xx = np.indices((512, 512)) - 256
rr = np.hypot(xx, yy)
xx = np.where(rr <= 256, xx, 0)
yy = np.where(rr <= 256, yy, 0)
zz = np.where(rr <= 256, 0, -1)
h, l, s = hls_from_vector(xx, yy, zz)
rgb = rgb_from_hls(h, l, s, mode=mode)
color_wheel = Image.fromarray(rgb)
# Plot the color wheel:
fig = plt.figure(figsize=(4, 4))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.imshow(color_wheel, origin='lower')
plt.tick_params(axis='both', which='both', labelleft='off', labelbottom='off',
left='off', right='off', top='off', bottom='off')
return
def rgb_from_hls(h, l, s, mode='triadic'):
"""Construct a rgb tuple from hue, luminance and saturation.
Parameters
----------
h : float, [0, 1)
Colorindex specifying the directional color according to the `mode` colormap scheme.
The colormap is periodic so that a value of 1 is equivalent to 0 again.
l : [0, 1] float
Luminance or lightness is the radiant power of the color. In the HLS model, it is defined
as the average between the largest and smallest color component. This definition puts
the primary (rgb) and secondary (ymc) colors into a plane halfway between black and
white for a luminance of 0.5. Directions pointing upwards will be encoded in white,
downwards in black.
s : [0, 1] float
The chroma scaled to fill the interval between 0 and 1. This means that smaller
amplitudes are encoded in grey, while large amplitudes will be in full color.
mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme of the color wheel. The default is the
standard HLS scheme with the three primary colors (red, blue, green) spaced with 120°
between them (triadic arrangement). The other option is a tetradic scheme commonly
used in holography where red, green, blue and yellow are used with 90° spacing. In both
cases, the saturation decreases to the center of the circle.
Returns
-------
rgb : tuple (N=3)
rgb tuple with the encoded direction color.
Notes
-----
Also works with numpy arrays as input! Always returns array of shape (..., 3)!
"""
_log.debug('Calling rgb_from_hls')
h = np.asarray(h)
l = np.asarray(l)
s = np.asarray(s)
hls_to_rgb = np.vectorize(colorsys.hls_to_rgb)
rgb_to_hls = np.vectorize(colorsys.rgb_to_hls)
if mode == 'triadic':
r, g, b = hls_to_rgb(h, l, s)
elif mode == 'tetradic':
rgb = hls_tetradic_cmap(h)
rh, gh, bh = rgb[..., 0], rgb[..., 1], rgb[..., 2]
rh = np.asarray(255 * rh).astype(np.uint8)
gh = np.asarray(255 * gh).astype(np.uint8)
bh = np.asarray(255 * bh).astype(np.uint8)
h_standard = rgb_to_hls(rh, gh, bh)[0]
r, g, b = hls_to_rgb(h_standard, l, s)
else:
raise ValueError('Mode not defined!')
return (255 * np.stack((r, g, b), axis=-1)).astype(np.uint8)
def hls_from_vector(x, y, z):
"""Construct a hls tuple from three coordinates representing a 3D direction.
Parameters
----------
x : float
x-coordinate of the desired direction to encode.
y : float
y-coordinate of the desired direction to encode.
z : float
z-coordinate of the desired direction to encode.
Returns
-------
hls : tuple (N=3)
hls tuple with the encoded direction color.
Notes
-----
Also works with numpy arrays as input!
"""
_log.debug('Calling hls_from_vector')
x = np.asarray(x)
y = np.asarray(y)
z = np.asarray(z)
phi = np.asarray(np.arctan2(y, x))
phi[phi < 0] += 2 * np.pi
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
theta = np.arccos(z / (r + 1E-30))
h = phi / (2 * np.pi)
l = 1 - theta / np.pi
s = r / r.max()
return h, l, s
def rgb_from_vector(x, y, z, mode='triadic'):
"""Construct a rgb tuple from three coordinates representing a 3D direction.
Parameters
----------
x : float
x-coordinate of the desired direction to encode.
y : float
y-coordinate of the desired direction to encode.
z : float
z-coordinate of the desired direction to encode.
mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme of the color wheel. The default is the
standard HLS scheme with the three primary colors (red, blue, green) spaced with 120°
between them (triadic arrangement). The other option is a tetradic scheme commonly
used in holography where red, green, blue and yellow are used with 90° spacing. In both
cases, the saturation decreases to the center of the circle.
Returns
-------
rgb : tuple (N=3)
rgb tuple with the encoded direction color.
Notes
-----
Also works with numpy arrays as input!
"""
_log.debug('Calling rgb_from_vector')
return rgb_from_hls(*hls_from_vector(x, y, z), mode=mode)
transparent_cmap = TransparentColormap(0.2, 0.3, 0.2, [0.75, 0.])
hls_tetradic_cmap = HLSTetradicColormap()
hls_triadic_cmap = HLSTriadicColormap()
pyramid/costfunction.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.Costfunction` class which represents a strategy to calculate
the so called `cost` of a threedimensional magnetization distribution."""
import logging
import numpy as np
from pyramid.regularisator import NoneRegularisator
__all__ = ['Costfunction']
class Costfunction(object):
"""Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
Represents a strategy for the calculation of the `cost` of a 3D magnetic distribution in
relation to two-dimensional phase maps. The `cost` is a measure for the difference of the
simulated phase maps from the magnetic distributions to the given set of phase maps and from
a priori knowledge represented by a :class:`~.Regularisator` object. Furthermore this class
provides convenient methods for the calculation of the derivative :func:`~.jac` or the product
with the Hessian matrix :func:`~.hess_dot` of the costfunction, which can be used by
optimizers. All required data should be given in a :class:`~DataSet` object.
Attributes
----------
fwd_model : :class:`~.ForwardModel`
The Forward model instance which should be used for the simulation of the phase maps which
will be compared to `y`.
regularisator : :class:`~.Regularisator`, optional
Regularisator class that's responsible for the regularisation term. If `None` or none is
given, no regularisation will be used.
y : :class:`~numpy.ndarray` (N=1)
Vector which lists all pixel values of all phase maps one after another.
m: int
Size of the image space.
n: int
Size of the input space.
Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `m x m` with m
being the length of the targetvector y.
"""
_log = logging.getLogger(__name__ + '.Costfunction')
def __init__(self, fwd_model, regularisator=None):
self._log.debug('Calling __init__')
self.fwd_model = fwd_model
if regularisator is None:
self.regularisator = NoneRegularisator()
else:
self.regularisator = regularisator
# Extract information from fwd_model:
self.y = self.fwd_model.y
self.n = self.fwd_model.n
self.m = self.fwd_model.m
self.Se_inv = self.fwd_model.Se_inv
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(fwd_model=%r, regularisator=%r)' % \
(self.__class__, self.fwd_model, self.regularisator)
def __str__(self):
self._log.debug('Calling __str__')
return 'Costfunction(fwd_model=%s, fwd_model=%s, regularisator=%s)' % \
(self.fwd_model, self.fwd_model, self.regularisator)
def __call__(self, x):
delta_y = self.fwd_model(x) - self.y
self.chisq_m = delta_y.dot(self.Se_inv.dot(delta_y))
self.chisq_a = self.regularisator(x)
self.chisq = self.chisq_m + self.chisq_a
return self.chisq
def init(self, x):
"""Initialise the costfunction by calculating the different cost terms.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution, for which the cost is calculated.
Returns
-------
None
"""
self._log.debug('Calling init')
self(x)
def jac(self, x):
"""Calculate the derivative of the costfunction for a given magnetization distribution.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution, for which the Jacobi vector is calculated.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Jacobi vector which represents the cost derivative of all voxels of the magnetization.
"""
assert len(x) == self.n
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x) - self.y))
+ self.regularisator.jac(x))
def hess_dot(self, x, vector):
"""Calculate the product of a `vector` with the Hessian matrix of the costfunction.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution at which the Hessian is calculated. The Hessian
is constant in this case, thus `x` can be set to None (it is not used int the
computation). It is implemented for the case that in the future nonlinear problems
have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution which is multiplied by the Hessian.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the input `vector` with the Hessian matrix of the costfunction.
"""
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
+ self.regularisator.hess_dot(x, vector))
def hess_diag(self, _):
""" Return the diagonal of the Hessian.
Parameters
----------
_ : undefined
Unused input
"""
return np.ones(self.n)
pyramid/dataset.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.DataSet` class for the collection of phase maps
and additional data like corresponding projectors."""
import logging
from numbers import Number
from pyramid.kernel import Kernel
from pyramid.phasemap import PhaseMap
from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.projector import Projector
import matplotlib.pyplot as plt
import numpy as np
from scipy import sparse
__all__ = ['DataSet']
class DataSet(object):
"""Class for collecting phase maps and corresponding projectors.
Represents a collection of (e.g. experimentally derived) phase maps, stored as
:class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
objects. At creation, the grid spacing `a` and the dimension `dim` of the magnetization
distribution have to be given. Data can be added via the :func:`~.append` method, where
a :class:`~.PhaseMap`, a :class:`~.Projector` and additional info have to be given.
Attributes
----------
a: float
The grid spacing in nm.
dim: tuple (N=3)
Dimensions of the 3D magnetization distribution.
b_0: double
The saturation induction in `T`.
mask: :class:`~numpy.ndarray` (N=3), optional
A boolean mask which defines the magnetized volume in 3D.
Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
being the length of the targetvector y (vectorized phase map information).
projectors: list of :class:`~.Projector`
A list of all stored :class:`~.Projector` objects.
phasemaps: list of :class:`~.PhaseMap`
A list of all stored :class:`~.PhaseMap` objects.
phase_vec: :class:`~numpy.ndarray` (N=1)
The concatenaded, vectorized phase of all :class:`~.PhaseMap` objects.
count(self): int
Number of phase maps and projectors in the dataset.
hook_points(self): :class:`~numpy.ndarray` (N=1)
Hook points which determine the start of values of a phase map in the `phase_vec`.
The length is `count + 1`.
"""
_log = logging.getLogger(__name__ + '.DataSet')
@property
def a(self):
"""The grid spacing in nm."""
return self._a
@a.setter
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = float(a)
@property
def mask(self):
"""A boolean mask which defines the magnetized volume in 3D."""
return self._mask
@mask.setter
def mask(self, mask):
if mask is not None:
assert mask.shape == self.dim, 'Mask dimensions must match!'
else:
mask = np.ones(self.dim, dtype=bool)
self._mask = mask.astype(np.bool)
@property
def m(self):
"""Size of the image space."""
return np.sum([len(p.phase_vec) for p in self.phasemaps])
@property
def n(self):
"""Size of the input space."""
return 3 * np.sum(self.mask)
@property
def count(self):
"""Number of phase maps and projectors in the dataset."""
return len(self.projectors)
@property
def phase_vec(self):
"""The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects."""
return np.concatenate([p.phase_vec for p in self.phasemaps])
@property
def hook_points(self):
"""Hook points which determine the start of values of a phase map in the `phase_vec`."""
result = [0]
for i, phasemap in enumerate(self.phasemaps):
result.append(result[i] + np.prod(phasemap.dim_uv))
return result
@property
def phasemappers(self):
"""List of phase mappers, created on demand with the projectors in mind."""
dim_uv_set = set([p.dim_uv for p in self.projectors])
kernel_list = [Kernel(self.a, dim_uv) for dim_uv in dim_uv_set]
return {kernel.dim_uv: PhaseMapperRDFC(kernel) for kernel in kernel_list}
def __init__(self, a, dim, b_0=1, mask=None, Se_inv=None):
self._log.debug('Calling __init__')
assert isinstance(dim, tuple) and len(dim) == 3, \
'Dimension has to be a tuple of length 3!'
self.a = a
self.dim = dim
self.b_0 = b_0
self.mask = mask
self.Se_inv = Se_inv
self.phasemaps = []
self.projectors = []
self._log.debug('Created: ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim=%r, b_0=%r, mask=%r, Se_inv=%r)' % (self.__class__, self.a, self.dim,
self.b_0, self.mask, self.Se_inv)
def __str__(self):
self._log.debug('Calling __str__')
return 'DataSet(a=%s, dim=%s, b_0=%s)' % (self.a, self.dim, self.b_0)
def append(self, phasemap, projector):
"""Appends a data pair of phase map and projection infos to the data collection.`
Parameters
----------
phasemap: :class:`~.PhaseMap`
A :class:`~.PhaseMap` object which should be added to the data collection.
projector: :class:`~.Projector`
A :class:`~.Projector` object which should be added to the data collection.
Returns
-------
None
"""
self._log.debug('Calling append')
assert isinstance(phasemap, PhaseMap) and isinstance(projector, Projector), \
'Argument has to be a tuple of a PhaseMap and a Projector object!'
assert projector.dim == self.dim, '3D dimensions must match!'
assert phasemap.dim_uv == projector.dim_uv, 'Projection dimensions (dim_uv) must match!'
self.phasemaps.append(phasemap)
self.projectors.append(projector)
def create_phasemaps(self, magdata):
"""Create a list of phasemaps with the projectors in the dataset for a given
:class:`~.VectorData` object.
Parameters
----------
magdata : :class:`~.VectorData`
Magnetic distribution to which the projectors of the dataset should be applied.
Returns
-------
phasemaps : list of :class:`~.phasemap.PhaseMap`
A list of the phase maps resulting from the projections specified in the dataset.
"""
self._log.debug('Calling create_phasemaps')
phasemaps = []
for projector in self.projectors:
mag_proj = projector(magdata)
phasemap = self.phasemappers[projector.dim_uv](mag_proj)
phasemap.mask = mag_proj.get_mask()[0, ...]
phasemaps.append(phasemap)
return phasemaps
def set_Se_inv_block_diag(self, cov_list):
"""Set the Se_inv matrix as a block diagonal matrix
Parameters
----------
cov_list: list of :class:`~numpy.ndarray`
List of inverted covariance matrices (one for each projection).
Returns
-------
None
"""
self._log.debug('Calling set_Se_inv_block_diag')
assert len(cov_list) == len(self.phasemaps), 'Needs one covariance matrix per phase map!'
self.Se_inv = sparse.block_diag(cov_list).tocsr()
def set_Se_inv_diag_with_conf(self, conf_list=None):
"""Set the Se_inv matrix as a block diagonal matrix from a list of confidence matrizes.
Parameters
----------
conf_list: list of :class:`~numpy.ndarray` (optional)
List of 2D confidence matrizes (one for each projection) which define trust regions.
If not given this uses the confidence matrizes of the phase maps.
Returns
-------
None
"""
self._log.debug('Calling set_Se_inv_diag_with_conf')
if conf_list is None: # if no confidence matrizes are given, extract from the phase maps!
conf_list = [phasemap.confidence for phasemap in self.phasemaps]
cov_list = [sparse.diags(c.ravel().astype(np.float32), 0) for c in conf_list]
self.set_Se_inv_block_diag(cov_list)
def set_3d_mask(self, mask_list=None):
"""Set the 3D mask from a list of 2D masks.
Parameters
----------
mask_list: list of :class:`~numpy.ndarray` (optional)
List of 2D masks, which represent the projections of the 3D mask. If not given this
uses the mask matrizes of the phase maps. If just one phase map is present, the
according mask is simply expanded to 3D and used directly.
Returns
-------
None
"""
self._log.debug('Calling set_3d_mask')
if mask_list is None: # if no masks are given, extract from phase maps:
mask_list = [phasemap.mask for phasemap in self.phasemaps]
if len(mask_list) == 1: # just one phasemap --> 3D mask equals 2D mask
self.mask = np.expand_dims(mask_list[0], axis=0) # z-dim is set to 1!
else: # 3D mask has to be constructed from 2D masks:
mask_3d_inv = np.zeros(self.dim)
for i, projector in enumerate(self.projectors):
mask_2d_inv = np.logical_not(self.phasemaps[i].mask.reshape(-1)) # inv. 2D mask
# Add extrusion of inv. 2D mask:
mask_3d_inv += projector.weight.T.dot(mask_2d_inv).reshape(self.dim)
self.mask = np.where(mask_3d_inv == 0, True, False)
def display_mask(self, ar_dens=1):
"""If it exists, display the 3D mask of the magnetization distribution.
Parameters
----------
ar_dens: int (optional)
Number defining the cell density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
Returns
-------
None
"""
self._log.debug('Calling display_mask')
if self.mask is not None:
from mayavi import mlab
zz, yy, xx = np.indices(self.dim)
ad = ar_dens
zz = zz[::ad, ::ad, ::ad].ravel()
yy = yy[::ad, ::ad, ::ad].ravel()
xx = xx[::ad, ::ad, ::ad].ravel()
mask_vec = self.mask[::ad, ::ad, ::ad].ravel().astype(dtype=np.int)
mlab.figure(size=(750, 700))
plot = mlab.points3d(xx, yy, zz, mask_vec, opacity=0.5,
mode='cube', scale_factor=ar_dens)
mlab.outline(plot)
mlab.axes(plot)
return plot
def phase_plots(self, magdata=None, title='Phase Map',
cmap='RdBu', limit=None, norm=None):
"""Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
Parameters
----------
magdata : :class:`~.VectorData`, optional
Magnetic distribution to which the projectors of the dataset should be applied. If not
given, the phasemaps in the dataset are used.
title : string, optional
The main part of the title of the plots. The default is 'Phase Map'. Additional
projector info is appended to this.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plots as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
Returns
-------
None
"""
self._log.debug('Calling phase_plots')
if magdata is not None:
phasemaps = self.create_phasemaps(magdata)
else:
phasemaps = self.phasemaps
[phasemap.plot_phase('{} ({})'.format(title, self.projectors[i].get_info()),
cmap=cmap, limit=limit, norm=norm)
for (i, phasemap) in enumerate(phasemaps)]
def combined_plots(self, magdata=None, title='Combined Plot', cmap='RdBu', limit=None,
norm=None, gain='auto', interpolation='none', grad_encode='bright'):
"""Display all phasemaps and the resulting color coded holography images.
Parameters
----------
magdata : :class:`~.VectorData`, optional
Magnetic distribution to which the projectors of the dataset should be applied. If not
given, the phasemaps in the dataset are used.
title : string, optional
The title of the plot. The default is 'Combined Plot'.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
gain : float, optional
The gain factor for determining the number of contour lines in the holographic
contour map. The default is 1.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
grad_encode: {'bright', 'dark', 'color', 'none'}, optional
Encoding mode of the phase gradient. 'none' produces a black-white image, 'color' just
encodes the direction (without gradient strength), 'dark' modulates the gradient
strength with a factor between 0 and 1 and 'bright' (which is the default) encodes
the gradient strength with color saturation.
Returns
-------
None
"""
self._log.debug('Calling combined_plots')
if magdata is not None:
phasemaps = self.create_phasemaps(magdata)
else:
phasemaps = self.phasemaps
for (i, phasemap) in enumerate(phasemaps):
phasemap.plot_combined('{} ({})'.format(title, self.projectors[i].get_info()),
cmap, limit, norm, gain, interpolation, grad_encode)
plt.show()
pyramid/diagnostics.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.Diagnostics` class for the calculation of diagnostics of a
specified costfunction for a fixed magnetization distribution."""
import logging
from pyramid import fft
from pyramid.fielddata import VectorData
from pyramid.phasemap import PhaseMap
import matplotlib.pyplot as plt
import numpy as np
import jutil
__all__ = ['Diagnostics']
class Diagnostics(object):
"""Class for calculating diagnostic properties of a specified costfunction.
For the calculation of diagnostic properties, a costfunction and a magnetization distribution
are specified at construction. With the :func:`~.set_position`, a position in 3D space can be
set at which all properties will be calculated. Properties are saved via boolean flags and
thus, calculation is only done if the position has changed in between. The standard deviation
and the measurement contribution require the execution of a conjugate gradient solver and can
take a while for larger problems.
Attributes
----------
x_rec: :class:`~numpy.ndarray`
Vectorized magnetization distribution at which the costfunction is evaluated.
cost: :class:`~.pyramid.costfunction.Costfunction`
Costfunction for which the diagnostics are calculated.
max_iter: int, optional
Maximum number of iterations. Default is 1000.
fwd_model: :class:`~pyramid.forwardmodel.ForwardModel`
Forward model used in the costfunction.
Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
being the length of the targetvector y (vectorized phase map information).
dim: tuple (N=3)
Dimensions of the 3D magnetic distribution.
row_idx: int
Row index of the system matrix corresponding to the current position in 3D space.
Notes
-----
Some properties depend on others, which may require recalculations of these prior
properties if necessary. The dependencies are ('-->' == 'requires'):
avrg_kern_row --> gain_row --> std --> m_inv_row
measure_contribution is independant
"""
_log = logging.getLogger(__name__ + '.Diagnostics')
@property
def cov_row(self):
"""Row of the covariance matrix (``S_a^-1+F'(x_f)^T S_e^-1 F'(x_f)``) which is needed for
the calculation of the gain and averaging kernel matrizes and which ideally contains the
variance at position `row_idx` for the current component and position in 3D."""
if not self._updated_cov_row:
e_i = fft.zeros(self.cost.n, dtype=fft.FLOAT)
e_i[self.row_idx] = 1
row = 2 * jutil.cg.conj_grad_solve(self._A, e_i, P=self._P, max_iter=self.max_iter)
self._std_row = np.asarray(row)
self._updated_cov_row = True
return self._std_row
@property
def std(self):
"""Standard deviation of the chosen component at the current position (calculated when
needed)."""
return np.sqrt(self.cov_row[self.row_idx])
@property
def gain_row(self):
"""Row of the gain matrix, which maps differences of phase measurements onto differences in
the retrieval result of the magnetization distribution(calculated when needed)."""
if not self._updated_gain_row:
self._gain_row = self.Se_inv.dot(self.fwd_model.jac_dot(self.x_rec, self.cov_row))
self._updated_gain_row = True
return self._gain_row
@property
def avrg_kern_row(self):
"""Row of the averaging kernel matrix (which is ideally the identity matrix), which
describes the smoothing introduced by the regularization (calculated when needed)."""
if not self._updated_avrg_kern_row:
self._avrg_kern_row = self.fwd_model.jac_T_dot(self.x_rec, self.gain_row)
self._updated_avrg_kern_row = True
return self._avrg_kern_row
@property
def measure_contribution(self):
"""The sum over an averaging kernel matrix row, which is an indicator for wheter a point of
the solution is determined by the measurement (close to `1`) or by a priori information
(close to `0`)."""
if not self._updated_measure_contribution:
cache = self.fwd_model.jac_dot(self.x_rec, fft.ones(self.cost.n, fft.FLOAT))
cache = self.fwd_model.jac_T_dot(self.x_rec, self.Se_inv.dot(cache))
mc = 2 * jutil.cg.conj_grad_solve(self._A, cache, P=self._P, max_iter=self.max_iter)
self._measure_contribution = mc
self._updated_measure_contribution = True
return self._measure_contribution
@property
def pos(self):
"""The current solution position, which specifies the 3D-point (and the component) of the
magnetization, for which diagnostics are calculated."""
return self._pos
@pos.setter
def pos(self, pos):
c, z, y, x = pos
assert self.mask[z, y, x], 'Position is outside of the provided mask!'
mask_vec = self.mask.ravel()
idx_3d = z * self.dim[1] * self.dim[2] + y * self.dim[2] + x
row_idx = c * np.prod(mask_vec.sum()) + mask_vec[:idx_3d].sum()
if row_idx != self.row_idx:
self._pos = pos
self.row_idx = row_idx
self._updated_cov_row = False
self._updated_gain_row = False
self._updated_avrg_kern_row = False
self._updated_measure_contribution = False
def __init__(self, x_rec, cost, max_iter=1000):
self._log.debug('Calling __init__')
self.x_rec = x_rec
self.cost = cost
self.max_iter = max_iter
self.fwd_model = cost.fwd_model
self.Se_inv = cost.Se_inv
self.dim = cost.data_set.dim
self.mask = cost.data_set.mask
self.row_idx = None
self.pos = (0,) + tuple(np.array(np.where(self.mask))[:, 0]) # first True mask entry
self._updated_cov_row = False
self._updated_gain_row = False
self._updated_avrg_kern_row = False
self._updated_measure_contribution = False
self._A = jutil.operator.CostFunctionOperator(self.cost, self.x_rec)
self._P = jutil.preconditioner.CostFunctionPreconditioner(self.cost, self.x_rec)
self._log.debug('Creating ' + str(self))
def get_avg_kern_row(self, pos=None):
"""Get the averaging kernel matrix row represented as a 3D magnetization distribution.
Parameters
----------
pos: tuple (N=4)
Position in 3D plus component `(c, z, y, x)`
Returns
-------
magdata_avg_kern: :class:`~pyramid.fielddata.VectorData`
Averaging kernel matrix row represented as a 3D magnetization distribution
"""
self._log.debug('Calling get_avg_kern_row')
if pos is not None:
self.pos = pos
magdata_avg_kern = VectorData(self.cost.data_set.a, fft.zeros((3,) + self.dim))
magdata_avg_kern.set_vector(self.avrg_kern_row, mask=self.mask)
return magdata_avg_kern
def calculate_averaging(self, pos=None):
"""Calculate and plot the averaging pixel number at a specified position for x, y or z.
Parameters
----------
pos: tuple (N=4)
Position in 3D plus component `(c, z, y, x)`
Returns
-------
px_avrg: float
The number of pixels over which is approximately averaged. The inverse is the FWHM.
(x, y, z): tuple of floats
The field of the averaging kernel summed along two axes (the remaining are x, y, z,
respectively).
Notes
-----
Uses the :func:`~.get_avg_kern_row` function
"""
self._log.debug('Calling calculate_averaging')
magdata_avg_kern = self.get_avg_kern_row(pos)
mag_x, mag_y, mag_z = magdata_avg_kern.field
x = mag_x.sum(axis=(0, 1))
y = mag_y.sum(axis=(0, 2))
z = mag_z.sum(axis=(1, 2))
# Plot helpful stuff:
plt.figure()
plt.axhline(y=0, ls='-', color='k')
plt.axhline(y=1, ls='-', color='k')
plt.plot(x, label='x', color='r', marker='o')
plt.plot(y, label='y', color='g', marker='o')
plt.plot(z, label='z', color='b', marker='o')
c = self.pos[0] # x, y or z component
data = [x, y, z][c]
col = ['r', 'g', 'b'][c]
i_m = np.argmax(data) # Index of the maximum
plt.axhline(y=data[i_m], ls='-', color=col)
plt.axhline(y=data[i_m] / 2, ls='--', color=col)
# Left side:
l = i_m
for i in np.arange(i_m - 1, -1, -1):
if data[i] < data[i_m] / 2:
# Linear interpolation between i and i + 1 to find left fractional index position:
l = (data[i_m] / 2 - data[i]) / (data[i + 1] - data[i]) + i
break
# Right side:
r = i_m
for i in np.arange(i_m + 1, data.size):
if data[i] < data[i_m] / 2:
# Linear interpolation between i and i - 1 to find right fractional index position:
r = (data[i_m] / 2 - data[i - 1]) / (data[i] - data[i - 1]) + i - 1
break
# Calculate FWHM:
fwhm = r - l
px_avrg = 1 / fwhm
plt.vlines(x=[l, r], ymin=0, ymax=data[i_m] / 2, linestyles=':', color=col)
# Add legend:
plt.legend()
return px_avrg, (x, y, z)
def get_gain_row_maps(self, pos=None):
"""Get the gain matrix row represented as a list of 2D (inverse) phase maps.
Parameters
----------
pos: tuple (N=4)
Position in 3D plus component `(c, z, y, x)`
Returns
-------
gain_map_list: list of :class:`~pyramid.phasemap.PhaseMap`
Gain matrix row represented as a list of 2D phase maps
Notes
-----
Note that the produced gain maps define the magnetization change at the current position
in 3d per phase change at the position of the . Take this into account when plotting the
maps (1/rad instead of rad).
"""
self._log.debug('Calling get_gain_row_maps')
if pos is not None:
self.pos = pos
hp = self.cost.data_set.hook_points
gain_map_list = []
for i, projector in enumerate(self.cost.data_set.projectors):
gain = self.gain_row[hp[i]:hp[i + 1]].reshape(projector.dim_uv)
gain_map_list.append(PhaseMap(self.cost.data_set.a, gain))
return gain_map_list
pyramid/fft.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Custom FFT module with numpy and FFTW support.
This module provides custom methods for FFTs including inverse, adjoint and real variants. The
FFTW library is supported and is used as a default if the import succeeds. Otherwise the numpy.fft
pack will be used. FFTW objects are saved in a cache after creation which speeds up further similar
FFT operations.
"""
import pickle
import logging
import os
import numpy as np
_log = logging.getLogger(__name__)
try:
import pyfftw
BACKEND = 'fftw'
except ImportError:
pyfftw = None
BACKEND = 'numpy'
_log.info('pyFFTW module not found. Using numpy implementation.')
try:
import multiprocessing
NTHREADS = multiprocessing.cpu_count()
del multiprocessing
except ImportError:
NTHREADS = 1
_log.info('multiprocessing module not found. Using single core.')
__all__ = ['plans', 'FLOAT', 'COMPLEX', 'dump_wisdom', 'load_wisdom',
'zeros', 'empty', 'ones', 'configure_backend',
'fftn', 'ifftn', 'rfftn', 'irfftn', 'rfftn_adj', 'irfftn_adj']
class FFTWCache(object):
"""Class for adding FFTW Plans and on-demand lookups.
This class is instantiated in this module to store FFTW plans and for the lookup of the former.
Attributes
----------
cache: dict
Cache for storing the FFTW plans.
Notes
-----
This class is used internally and is not normally not intended to be used directly by the user.
"""
_log = logging.getLogger(__name__ + '.FFTWCache')
def __init__(self):
self._log.debug('Calling __init__')
self.cache = dict()
self._log.debug('Created ' + str(self))
def add_fftw(self, fft_type, fftw_obj, s, axes, nthreads):
"""Add an FFTW object to the cache.
Parameters
----------
fft_type: basestring
Identifier sting for the FFT type ('fftn', 'ifftn', 'rfftn', 'irfftn').
fftw_obj: :class:`~pyfftw.FFTW` object
The FFTW object which should be added to the cache.
s: tuple of ints
Shape of the output array.
axes: tuple of ints
The axes along which the FFTW should be executed.
nthreads: int
Number of threads which should be used.
"""
self._log.debug('Calling add_fftw')
in_arr = fftw_obj.get_input_array()
key = (fft_type, in_arr.shape, in_arr.dtype, s, axes, nthreads)
self.cache[key] = fftw_obj
def lookup_fftw(self, fft_type, in_arr, s, axes, nthreads):
"""
Parameters
----------
fft_type: basestring
Identifier sting for the FFT type ('fftn', 'ifftn', 'rfftn', 'irfftn').
in_arr:
Input array, internally, just the `dtype` and the `shape` are used to identify the FFT.
s: tuple of ints
Shape of the output array.
axes: tuple of ints
The axes along which the FFTW should be executed.
nthreads: int
Number of threads which should be used.
Returns
-------
fftw_obj: :class:`~pyfftw.FFTW` object
The requested FFTW object.
"""
self._log.debug('Calling lookup_fftw')
key = (fft_type, in_arr.shape, in_arr.dtype, s, axes, nthreads)
return self.cache.get(key, None)
def clear_cache(self):
"""Clear the cache."""
self._log.debug('Calling clear_cache')
self.cache = dict()
plans = FFTWCache()
FLOAT = np.float32 # One convenient place to
COMPLEX = np.complex64 # change from 32 to 64 bit
# Numpy functions:
def _fftn_numpy(a, s=None, axes=None):
return np.fft.fftn(a, s, axes)
def _ifftn_numpy(a, s=None, axes=None):
return np.fft.ifftn(a, s, axes)
def _rfftn_numpy(a, s=None, axes=None):
return np.fft.rfftn(a, s, axes)
def _irfftn_numpy(a, s=None, axes=None):
return np.fft.irfftn(a, s, axes)
def _rfftn_adj_numpy(a):
n = 2 * (a.shape[-1] - 1)
out_shape = a.shape[:-1] + (n,)
out_arr = zeros(out_shape, dtype=a.dtype)
out_arr[:, :n] = a
return _ifftn_numpy(out_arr).real * np.prod(out_shape)
def _irfftn_adj_numpy(a):
n = a.shape[-1] // 2 + 1
out_arr = _fftn_numpy(a, axes=(-1,)) / a.shape[-1]
if a.shape[-1] % 2 == 0: # even
out_arr[:, 1:n - 1] += np.conj(out_arr[:, :n - 1:-1])
else: # odd
out_arr[:, 1:n] += np.conj(out_arr[:, :n - 1:-1])
axes = tuple(range(len(out_arr.shape[:-1])))
return _fftn_numpy(out_arr[:, :n], axes=axes) / np.prod(out_arr.shape[:-1])
# FFTW functions:
def _fftn_fftw(a, s=None, axes=None):
fftw = plans.lookup_fftw('fftn', a, s, axes, NTHREADS)
if fftw is None:
fftw = pyfftw.builders.fftn(a, s, axes, threads=NTHREADS)
plans.add_fftw('fftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
def _ifftn_fftw(a, s=None, axes=None):
fftw = plans.lookup_fftw('ifftn', a, s, axes, NTHREADS)
if fftw is None:
fftw = pyfftw.builders.ifftn(a, s, axes, threads=NTHREADS)
plans.add_fftw('ifftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
def _rfftn_fftw(a, s=None, axes=None):
fftw = plans.lookup_fftw('rfftn', a, s, axes, NTHREADS)
if fftw is None:
fftw = pyfftw.builders.rfftn(a, s, axes, threads=NTHREADS)
plans.add_fftw('rfftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
def _irfftn_fftw(a, s=None, axes=None):
fftw = plans.lookup_fftw('irfftn', a, s, axes, NTHREADS)
if fftw is None:
fftw = pyfftw.builders.irfftn(a, s, axes, threads=NTHREADS)
plans.add_fftw('irfftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
def _rfftn_adj_fftw(a):
# Careful: just works for even a (which is guaranteed by the kernel!)
n = 2 * (a.shape[-1] - 1)
out_shape = a.shape[:-1] + (n,)
out_arr = zeros(out_shape, dtype=a.dtype)
out_arr[:, :a.shape[-1]] = a
return _ifftn_fftw(out_arr).real * np.prod(out_shape)
def _irfftn_adj_fftw(a):
out_arr = _fftn_fftw(a, axes=(-1,)) / a.shape[-1] # FFT of last axis
n = a.shape[-1] // 2 + 1
if a.shape[-1] % 2 == 0: # even
out_arr[:, 1:n - 1] += np.conj(out_arr[:, :n - 1:-1])
else: # odd
out_arr[:, 1:n] += np.conj(out_arr[:, :n - 1:-1])
axes = tuple(range(len(out_arr.shape[:-1])))
return _fftn_fftw(out_arr[:, :n], axes=axes) / np.prod(out_arr.shape[:-1])
# These wisdom functions do nothing if pyFFTW is not available:
def dump_wisdom(fname):
"""Wrapper function for the pyfftw.export_wisdom(), which uses a pickle dump.
Parameters
----------
fname: string
Name of the file in which the wisdom is saved.
Returns
-------
None
"""
_log.debug('Calling dump_wisdom')
if pyfftw is not None:
with open(fname, 'wb') as fp:
pickle.dump(pyfftw.export_wisdom(), fp, pickle.HIGHEST_PROTOCOL)
def load_wisdom(fname):
"""Wrapper function for the pyfftw.import_wisdom(), which uses a pickle to load a file.
Parameters
----------
fname: string
Name of the file from which the wisdom is loaded.
Returns
-------
None
"""
_log.debug('Calling load_wisdom')
if pyfftw is not None:
if not os.path.exists(fname):
print("Warning: Wisdom file does not exist. First time use?")
else:
with open(fname, 'rb') as fp:
pyfftw.import_wisdom(pickle.load(fp))
# Array setups:
def empty(shape, dtype=FLOAT):
"""Return a new array of given shape and type without initializing entries.
Parameters
----------
shape: int or tuple of int
Shape of the array.
dtype: data-type, optional
Desired output data-type.
Returns
-------
out: :class:`~numpy.ndarray`
The created array.
"""
_log.debug('Calling empty')
result = np.empty(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
return result
def zeros(shape, dtype=FLOAT):
"""Return a new array of given shape and type, filled with zeros.
Parameters
----------
shape: int or tuple of int
Shape of the array.
dtype: data-type, optional
Desired output data-type.
Returns
-------
out: :class:`~numpy.ndarray`
The created array.
"""
_log.debug('Calling zeros')
result = np.zeros(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
return result
def ones(shape, dtype=FLOAT):
"""Return a new array of given shape and type, filled with ones.
Parameters
----------
shape: int or tuple of int
Shape of the array.
dtype: data-type, optional
Desired output data-type.
Returns
-------
out: :class:`~numpy.ndarray`
The created array.
"""
_log.debug('Calling ones')
result = np.ones(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
return result
# Configure backend:
def configure_backend(backend):
"""Change FFT backend.
Parameters
----------
backend: string
Backend to use. Supported values are "numpy" and "fftw".
Returns
-------
None
"""
_log.debug('Calling configure_backend')
global fftn
global ifftn
global rfftn
global irfftn
global rfftn_adj
global irfftn_adj
global BACKEND
if backend == 'numpy':
fftn = _fftn_numpy
ifftn = _ifftn_numpy
rfftn = _rfftn_numpy
irfftn = _irfftn_numpy
rfftn_adj = _rfftn_adj_numpy
irfftn_adj = _irfftn_adj_numpy
BACKEND = 'numpy'
elif backend == 'fftw':
if pyfftw is not None:
fftn = _fftn_fftw
ifftn = _ifftn_fftw
rfftn = _rfftn_fftw
irfftn = _irfftn_fftw
rfftn_adj = _rfftn_adj_fftw
irfftn_adj = _irfftn_adj_fftw
BACKEND = 'pyfftw'
else:
print('Error: FFTW requested but not available')
# On import:
ifftn = None
fftn = None
rfftn = None
irfftn = None
rfftn_adj = None
irfftn_adj = None
if pyfftw is not None:
configure_backend('fftw')
else:
configure_backend('numpy')
pyramid/fieldconverter.py deleted 100644 → 0
View file @ 68453aaa
# coding=utf-8
"""Convert vector fields.
The :mod:`~.fieldconverter` provides methods for converting a magnetization distribution `M` into
a vector potential `A` and convert this in turn into a magnetic field `B`. The direct way is also
possible.
"""
import logging
import numpy as np
from pyramid import fft
from pyramid.fielddata import VectorData
__all__ = ['convert_M_to_A', 'convert_A_to_B', 'convert_M_to_B']
_log = logging.getLogger(__name__)
def convert_M_to_A(magdata, b_0=1.0):
"""Convert a magnetic vector distribution into a vector potential `A`.
Parameters
----------
magdata: :class:`~pyramid.magdata.VectorData` object
The magnetic vector field from which the A-field is calculated.
b_0: float, optional
The saturation magnetization which is used in the calculation.
Returns
-------
b_data: :class:`~pyramid.magdata.VectorData` object
The calculated B-field.
"""
_log.debug('Calling convert_M_to_A')
# Preparations of variables:
assert isinstance(magdata, VectorData), 'Only VectorData objects can be mapped!'
dim = magdata.dim
dim_kern = tuple(2 * np.array(dim) - 1) # Dimensions of the kernel
if fft.BACKEND == 'pyfftw':
dim_pad = tuple(2 * np.array(dim)) # is at least even (not neccessary a power of 2)
elif fft.BACKEND == 'numpy':
dim_pad = tuple(2 ** np.ceil(np.log2(2 * np.array(dim))).astype(int)) # pow(2)
else:
raise ValueError('Backend of the fft module is not correctly initiated!')
slice_B = (slice(dim[0] - 1, dim_kern[0]), # Shift because kernel center
slice(dim[1] - 1, dim_kern[1]), # is not at (0, 0, 0)!
slice(dim[2] - 1, dim_kern[2]))
slice_M = (slice(0, dim[0]), # Magnetization is padded on the far end!
slice(0, dim[1]), # B-field cutout is shifted as listed above
slice(0, dim[2])) # because of the kernel center!
# Set up kernels
coeff = magdata.a * b_0 / (4 * np.pi)
zzz, yyy, xxx = np.indices(dim_kern)
xxx -= dim[2] - 1
yyy -= dim[1] - 1
zzz -= dim[0] - 1
k_x = fft.empty(dim_kern, dtype=fft.FLOAT)
k_y = fft.empty(dim_kern, dtype=fft.FLOAT)
k_z = fft.empty(dim_kern, dtype=fft.FLOAT)
k_x[...] = coeff * xxx / np.abs(xxx ** 2 + yyy ** 2 + zzz ** 2 + 1E-30) ** 3
k_y[...] = coeff * yyy / np.abs(xxx ** 2 + yyy ** 2 + zzz ** 2 + 1E-30) ** 3
k_z[...] = coeff * zzz / np.abs(xxx ** 2 + yyy ** 2 + zzz ** 2 + 1E-30) ** 3
# Calculate Fourier trafo of kernel components:
k_x_fft = fft.rfftn(k_x, dim_pad)
k_y_fft = fft.rfftn(k_y, dim_pad)
k_z_fft = fft.rfftn(k_z, dim_pad)
# Prepare magnetization:
x_mag = fft.zeros(dim_pad, dtype=fft.FLOAT)
y_mag = fft.zeros(dim_pad, dtype=fft.FLOAT)
z_mag = fft.zeros(dim_pad, dtype=fft.FLOAT)
x_mag[slice_M] = magdata.field[0, ...]
y_mag[slice_M] = magdata.field[1, ...]
z_mag[slice_M] = magdata.field[2, ...]
# Calculate Fourier trafo of magnetization components:
x_mag_fft = fft.rfftn(x_mag)
y_mag_fft = fft.rfftn(y_mag)
z_mag_fft = fft.rfftn(z_mag)
# Convolve:
a_x_fft = y_mag_fft * k_z_fft - z_mag_fft * k_y_fft
a_y_fft = z_mag_fft * k_x_fft - x_mag_fft * k_z_fft
a_z_fft = x_mag_fft * k_y_fft - y_mag_fft * k_x_fft
a_x = fft.irfftn(a_x_fft)[slice_B]
a_y = fft.irfftn(a_y_fft)[slice_B]
a_z = fft.irfftn(a_z_fft)[slice_B]
# Return A-field:
return VectorData(magdata.a, np.asarray((a_x, a_y, a_z)))
def convert_A_to_B(a_data):
"""Convert a vector potential `A` into a B-field distribution.
Parameters
----------
a_data: :class:`~pyramid.magdata.VectorData` object
The vector potential field from which the A-field is calculated.
Returns
-------
b_data: :class:`~pyramid.magdata.VectorData` object
The calculated B-field.
"""
_log.debug('Calling convert_A_to_B')
assert isinstance(a_data, VectorData), 'Only VectorData objects can be mapped!'
#
axis = tuple([i for i in range(3) if a_data.dim[i] > 1])
#
x_grads = np.gradient(a_data.field[0, ...], axis=axis) #/ a_data.a
y_grads = np.gradient(a_data.field[1, ...], axis=axis) #/ a_data.a
z_grads = np.gradient(a_data.field[2, ...], axis=axis) #/ a_data.a
#
x_gradii = np.zeros(a_data.shape)
y_gradii = np.zeros(a_data.shape)
z_gradii = np.zeros(a_data.shape)
#
for i, axis in enumerate(axis):
x_gradii[axis] = x_grads[i]
y_gradii[axis] = y_grads[i]
z_gradii[axis] = z_grads[i]
#
x_grad_z, x_grad_y, x_grad_x = x_gradii
y_grad_z, y_grad_y, y_grad_x = y_gradii
z_grad_z, z_grad_y, z_grad_x = z_gradii
# Calculate cross product:
b_x = (z_grad_y - y_grad_z)
b_y = (x_grad_z - z_grad_x)
b_z = (y_grad_x - x_grad_y)
# Return B-field:
return VectorData(a_data.a, np.asarray((b_x, b_y, b_z)))
# Calculate gradients:
x_mag, y_mag, z_mag = a_data.field
x_grad_z, x_grad_y, x_grad_x = np.gradient(x_mag)
y_grad_z, y_grad_y, y_grad_x = np.gradient(y_mag)
z_grad_z, z_grad_y, z_grad_x = np.gradient(z_mag)
# Calculate cross product:
b_x = (z_grad_y - y_grad_z)
b_y = (x_grad_z - z_grad_x)
b_z = (y_grad_x - x_grad_y)
# Return B-field:
return VectorData(a_data.a, np.asarray((b_x, b_y, b_z)))
def convert_M_to_B(magdata, b_0=1.0):
"""Convert a magnetic vector distribution into a B-field distribution.
Parameters
----------
magdata: :class:`~pyramid.magdata.VectorData` object
The magnetic vector field from which the B-field is calculated.
b_0: float, optional
The saturation magnetization which is used in the calculation.
Returns
-------
b_data: :class:`~pyramid.magdata.VectorData` object
The calculated B-field.
"""
_log.debug('Calling convert_M_to_B')
assert isinstance(magdata, VectorData), 'Only VectorData objects can be mapped!'
return convert_A_to_B(convert_M_to_A(magdata, b_0=b_0))
pyramid/fielddata.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides classes for storing vector and scalar 3D-field."""
import logging
import abc
from numbers import Number
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.ticker import MaxNLocator
from PIL import Image
from scipy.ndimage.interpolation import zoom
from . import fft
from . import colors
__all__ = ['VectorData', 'ScalarData']
class FieldData(object, metaclass=abc.ABCMeta):
"""Class for storing field data.
Abstract base class for the representatio of magnetic or electric fields (see subclasses).
Fields can be accessed as 3D numpy arrays via the `field` property or as a vector via
`field_vec`. :class:`~.FieldData` objects support negation, arithmetic operators
(``+``, ``-``, ``*``) and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers
and other :class:`~.FieldData` objects of the same subclass, if their dimensions and grid
spacings match. It is possible to load data from HDF5 or LLG (.txt) files or to save the data
in these formats. Specialised plotting methods are also provided.
Attributes
----------
a: float
The grid spacing in nm.
field: :class:`~numpy.ndarray` (N=4)
The field distribution for every 3D-gridpoint.
"""
_log = logging.getLogger(__name__ + '.FieldData')
@property
def a(self):
"""The grid spacing in nm."""
return self._a
@a.setter
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = float(a)
@property
def shape(self):
"""The shape of the `field` (3D for scalar, 4D vor vector field)."""
return self.field.shape
@property
def dim(self):
"""Dimensions (z, y, x) of the grid, only 3D coordinates, without components if present."""
return self.shape[-3:]
@property
def field(self):
"""The field strength for every 3D-gridpoint (scalar: 3D, vector: 4D)."""
return self._field
@field.setter
def field(self, field):
assert isinstance(field, np.ndarray), 'Field has to be a numpy array!'
assert 3 <= len(field.shape) <= 4, 'Field has to be 3- or 4-dimensional (scalar / vector)!'
if len(field.shape) == 4:
assert field.shape[0] == 3, 'A vector field has to have exactly 3 components!'
self._field = np.asarray(field, dtype=fft.FLOAT)
@property
def field_amp(self):
"""The field amplitude (returns the field itself for scalar and the vector amplitude
calculated via a square sum for a vector field."""
if len(self.shape) == 4:
return np.sqrt(np.sum(self.field ** 2, axis=0))
else:
return self.field
@property
def field_vec(self):
"""Vector containing the vector field distribution."""
return np.reshape(self.field, -1)
@field_vec.setter
def field_vec(self, mag_vec):
mag_vec = np.asarray(mag_vec, dtype=fft.FLOAT)
assert np.size(mag_vec) == np.prod(self.shape), \
'Vector has to match field shape! {} {}'.format(mag_vec.shape, np.prod(self.shape))
self.field = mag_vec.reshape((3,) + self.dim)
def __init__(self, a, field):
self._log.debug('Calling __init__')
self.a = a
self.field = field
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, field=%r)' % (self.__class__, self.a, self.field)
def __str__(self):
self._log.debug('Calling __str__')
return '%s(a=%s, dim=%s)' % (self.__class__, self.a, self.dim)
def __neg__(self): # -self
self._log.debug('Calling __neg__')
return self.__class__(self.a, -self.field)
def __add__(self, other): # self + other
self._log.debug('Calling __add__')
assert isinstance(other, (FieldData, Number)), \
'Only FieldData objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, Number): # other is a Number
self._log.debug('Adding an offset')
return self.__class__(self.a, self.field + other)
elif isinstance(other, FieldData):
self._log.debug('Adding two FieldData objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.shape == self.shape, 'Added field has to have the same dimensions!'
return self.__class__(self.a, self.field + other.field)
def __sub__(self, other): # self - other
self._log.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
self._log.debug('Calling __mul__')
assert isinstance(other, Number), 'FieldData objects can only be multiplied by numbers!'
return self.__class__(self.a, self.field * other)
def __truediv__(self, other): # self / other
self._log.debug('Calling __truediv__')
assert isinstance(other, Number), 'FieldData objects can only be divided by numbers!'
return self.__class__(self.a, self.field / other)
def __floordiv__(self, other): # self // other
self._log.debug('Calling __floordiv__')
assert isinstance(other, Number), 'FieldData objects can only be divided by numbers!'
return self.__class__(self.a, self.field // other)
def __radd__(self, other): # other + self
self._log.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
self._log.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
self._log.debug('Calling __rmul__')
return self.__mul__(other)
def __iadd__(self, other): # self += other
self._log.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
self._log.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
self._log.debug('Calling __imul__')
return self.__mul__(other)
def __itruediv__(self, other): # self /= other
self._log.debug('Calling __itruediv__')
return self.__truediv__(other)
def __ifloordiv__(self, other): # self //= other
self._log.debug('Calling __ifloordiv__')
return self.__floordiv__(other)
def __array__(self, dtype=None):
if dtype:
return self.field.astype(dtype)
else:
return self.field
def __array_wrap__(self, array, _=None): # _ catches the context, which is not used.
return type(self)(self.a, array)
def copy(self):
"""Returns a copy of the :class:`~.FieldData` object
Returns
-------
field_data: :class:`~.FieldData`
A copy of the :class:`~.FieldData`.
"""
self._log.debug('Calling copy')
return self.__class__(self.a, self.field.copy())
def get_mask(self, threshold=0):
"""Mask all pixels where the amplitude of the field lies above `threshold`.
Parameters
----------
threshold : float, optional
A pixel only gets masked, if it lies above this threshold . The default is 0.
Returns
-------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Mask of the pixels where the amplitude of the field lies above `threshold`.
"""
self._log.debug('Calling get_mask')
return np.where(self.field_amp > threshold, True, False)
def contour_plot3d(self, title='Field Distribution', contours=10, opacity=0.25):
"""Plot the field as a 3D-contour plot.
Parameters
----------
title: string, optional
The title for the plot.
contours: int, optional
Number of contours which should be plotted.
opacity: float, optional
Defines the opacity of the contours. Default is 0.25.
Returns
-------
plot : :class:`mayavi.modules.vectors.Vectors`
The plot object.
"""
self._log.debug('Calling quiver_plot3D')
from mayavi import mlab
# Plot them as vectors:
mlab.figure(size=(750, 700))
plot = mlab.contour3d(self.field_amp, contours=contours, opacity=opacity)
mlab.outline(plot)
mlab.axes(plot)
mlab.title(title, height=0.95, size=0.35)
mlab.orientation_axes()
return plot
@abc.abstractmethod
def scale_down(self, n):
"""Scale down the field distribution by averaging over two pixels along each axis.
Parameters
----------
n : int, optional
Number of times the field distribution is scaled down. The default is 1.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
Only possible, if each axis length is a power of 2!
"""
pass
@abc.abstractmethod
def scale_up(self, n, order):
"""Scale up the field distribution using spline interpolation of the requested order.
Parameters
----------
n : int, optional
Power of 2 with which the grid is scaled. Default is 1, which means every axis is
increased by a factor of ``2**1 = 2``.
order : int, optional
The order of the spline interpolation, which has to be in the range between 0 and 5
and defaults to 0.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
"""
pass
@abc.abstractmethod
def get_vector(self, mask):
"""Returns the field as a vector, specified by a mask.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Masks the pixels from which the entries should be taken.
Returns
-------
vector : :class:`~numpy.ndarray` (N=1)
The vector containing the field of the specified pixels.
"""
pass
@abc.abstractmethod
def set_vector(self, vector, mask):
"""Set the field of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean), optional
Masks the pixels from which the field should be taken.
vector : :class:`~numpy.ndarray` (N=1)
The vector containing the field of the specified pixels.
Returns
-------
None
"""
pass
@classmethod
def from_signal(cls, signal):
"""Convert a :class:`~hyperspy.signals.Signal` object to a :class:`~.FieldData` object.
Parameters
----------
signal: :class:`~hyperspy.signals.Signal`
The :class:`~hyperspy.signals.Signal` object which should be converted to FieldData.
Returns
-------
magdata: :class:`~.FieldData`
A :class:`~.FieldData` object containing the loaded data.
Notes
-----
This method recquires the hyperspy package!
"""
cls._log.debug('Calling from_signal')
return cls(signal.axes_manager[0].scale, signal.data)
@abc.abstractmethod
def to_signal(self):
"""Convert :class:`~.FieldData` data into a HyperSpy signal.
Returns
-------
signal: :class:`~hyperspy.signals.Signal`
Representation of the :class:`~.FieldData` object as a HyperSpy Signal.
Notes
-----
This method recquires the hyperspy package!
"""
self._log.debug('Calling to_signal')
try: # Try importing HyperSpy:
import hyperspy.api as hs
except ImportError:
self._log.error('This method recquires the hyperspy package!')
return
# Create signal:
signal = hs.signals.BaseSignal(self.field) # All axes are signal axes!
# Set axes:
signal.axes_manager[0].name = 'x-axis'
signal.axes_manager[0].units = 'nm'
signal.axes_manager[0].scale = self.a
signal.axes_manager[1].name = 'y-axis'
signal.axes_manager[1].units = 'nm'
signal.axes_manager[1].scale = self.a
signal.axes_manager[2].name = 'z-axis'
signal.axes_manager[2].units = 'nm'
signal.axes_manager[2].scale = self.a
return signal
class VectorData(FieldData):
"""Class for storing vector ield data.
Represents 3-dimensional vector field distributions with 3 components which are stored as a
3-dimensional numpy array in `field`, but which can also be accessed as a vector via
`field_vec`. :class:`~.VectorData` objects support negation, arithmetic operators
(``+``, ``-``, ``*``) and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers
and other :class:`~.VectorData` objects, if their dimensions and grid spacings match. It is
possible to load data from HDF5 or LLG (.txt) files or to save the data in these formats.
Plotting methods are also provided.
Attributes
----------
a: float
The grid spacing in nm.
field: :class:`~numpy.ndarray` (N=4)
The `x`-, `y`- and `z`-component of the vector field for every 3D-gridpoint
as a 4-dimensional numpy array (first dimension has to be 3, because of the 3 components).
"""
_log = logging.getLogger(__name__ + '.VectorData')
def scale_down(self, n=1):
"""Scale down the field distribution by averaging over two pixels along each axis.
Parameters
----------
n : int, optional
Number of times the field distribution is scaled down. The default is 1.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
Only possible, if each axis length is a power of 2!
"""
self._log.debug('Calling scale_down')
assert n > 0 and isinstance(n, int), 'n must be a positive integer!'
self.a *= 2 ** n
for t in range(n):
# Pad if necessary:
pz, py, px = self.dim[0] % 2, self.dim[1] % 2, self.dim[2] % 2
if pz != 0 or py != 0 or px != 0:
self.field = np.pad(self.field, ((0, 0), (0, pz), (0, py), (0, px)),
mode='constant')
# Create coarser grid for the vector field:
shape_4d = (3, self.dim[0] // 2, 2, self.dim[1] // 2, 2, self.dim[2] // 2, 2)
self.field = self.field.reshape(shape_4d).mean(axis=(6, 4, 2))
def scale_up(self, n=1, order=0):
"""Scale up the field distribution using spline interpolation of the requested order.
Parameters
----------
n : int, optional
Power of 2 with which the grid is scaled. Default is 1, which means every axis is
increased by a factor of ``2**1 = 2``.
order : int, optional
The order of the spline interpolation, which has to be in the range between 0 and 5
and defaults to 0.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
"""
self._log.debug('Calling scale_up')
assert n > 0 and isinstance(n, int), 'n must be a positive integer!'
assert 5 > order >= 0 and isinstance(order, int), \
'order must be a positive integer between 0 and 5!'
self.a /= 2 ** n
self.field = np.array((zoom(self.field[0], zoom=2 ** n, order=order),
zoom(self.field[1], zoom=2 ** n, order=order),
zoom(self.field[2], zoom=2 ** n, order=order)))
def pad(self, pad_values):
"""Pad the current field distribution with zeros for each individual axis.
Parameters
----------
pad_values : tuple of int
Number of zeros which should be padded. Provided as a tuple where each entry
corresponds to an axis. An entry can be one int (same padding for both sides) or again
a tuple which specifies the pad values for both sides of the corresponding axis.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions accordingly.
"""
self._log.debug('Calling pad')
assert len(pad_values) == 3, 'Pad values for each dimension have to be provided!'
pv = np.zeros(6, dtype=np.int)
for i, values in enumerate(pad_values):
assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
pv[2 * i:2 * (i + 1)] = values
self.field = np.pad(self.field, ((0, 0), (pv[0], pv[1]), (pv[2], pv[3]), (pv[4], pv[5])),
mode='constant')
def crop(self, crop_values):
"""Crop the current field distribution with zeros for each individual axis.
Parameters
----------
crop_values : tuple of int
Number of zeros which should be cropped. Provided as a tuple where each entry
corresponds to an axis. An entry can be one int (same cropping for both sides) or again
a tuple which specifies the crop values for both sides of the corresponding axis.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions accordingly.
"""
self._log.debug('Calling crop')
assert len(crop_values) == 3, 'Crop values for each dimension have to be provided!'
cv = np.zeros(6, dtype=np.int)
for i, values in enumerate(crop_values):
assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
cv[2 * i:2 * (i + 1)] = values
cv *= np.resize([1, -1], len(cv))
cv = np.where(cv == 0, None, cv)
self.field = self.field[:, cv[0]:cv[1], cv[2]:cv[3], cv[4]:cv[5]]
def get_vector(self, mask):
"""Returns the vector field components arranged in a vector, specified by a mask.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Masks the pixels from which the components should be taken.
Returns
-------
vector : :class:`~numpy.ndarray` (N=1)
The vector containing vector field components of the specified pixels.
Order is: first all `x`-, then all `y`-, then all `z`-components.
"""
self._log.debug('Calling get_vector')
if mask is not None:
return np.reshape([self.field[0][mask],
self.field[1][mask],
self.field[2][mask]], -1)
else:
return self.field_vec
def set_vector(self, vector, mask=None):
"""Set the field components of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean), optional
Masks the pixels from which the components should be taken.
vector : :class:`~numpy.ndarray` (N=1)
The vector containing vector field components of the specified pixels.
Order is: first all `x`-, then all `y-, then all `z`-components.
Returns
-------
None
"""
self._log.debug('Calling set_vector')
vector = np.asarray(vector, dtype=fft.FLOAT)
assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
count = np.size(vector) // 3
if mask is not None:
self.field[0][mask] = vector[:count] # x-component
self.field[1][mask] = vector[count:2 * count] # y-component
self.field[2][mask] = vector[2 * count:] # z-component
else:
self.field_vec = vector
def flip(self, axis='x'):
"""Flip/mirror the vector field around the specified axis.
Parameters
----------
axis: {'x', 'y', 'z'}, optional
The axis around which the vector field is flipped.
Returns
-------
magdata_flip: :class:`~.VectorData`
A flipped copy of the :class:`~.VectorData` object.
"""
self._log.debug('Calling flip')
if axis == 'x':
mag_x, mag_y, mag_z = self.field[:, :, :, ::-1]
field_flip = np.array((-mag_x, mag_y, mag_z))
elif axis == 'y':
mag_x, mag_y, mag_z = self.field[:, :, ::-1, :]
field_flip = np.array((mag_x, -mag_y, mag_z))
elif axis == 'z':
mag_x, mag_y, mag_z = self.field[:, ::-1, :, :]
field_flip = np.array((mag_x, mag_y, -mag_z))
else:
raise ValueError("Wrong input! 'x', 'y', 'z' allowed!")
return VectorData(self.a, field_flip)
def rot90(self, axis='x'):
"""Rotate the vector field 90° around the specified axis (right hand rotation).
Parameters
----------
axis: {'x', 'y', 'z'}, optional
The axis around which the vector field is rotated.
Returns
-------
magdata_rot: :class:`~.VectorData`
A rotated copy of the :class:`~.VectorData` object.
"""
self._log.debug('Calling rot90')
if axis == 'x':
field_rot = np.zeros((3, self.dim[1], self.dim[0], self.dim[2]))
for i in range(self.dim[2]):
mag_x, mag_y, mag_z = self.field[:, :, :, i]
mag_xrot, mag_yrot, mag_zrot = np.rot90(mag_x), np.rot90(mag_y), np.rot90(mag_z)
field_rot[:, :, :, i] = np.array((mag_xrot, mag_zrot, -mag_yrot))
elif axis == 'y':
field_rot = np.zeros((3, self.dim[2], self.dim[1], self.dim[0]))
for i in range(self.dim[1]):
mag_x, mag_y, mag_z = self.field[:, :, i, :]
mag_xrot, mag_yrot, mag_zrot = np.rot90(mag_x), np.rot90(mag_y), np.rot90(mag_z)
field_rot[:, :, i, :] = np.array((mag_zrot, mag_yrot, -mag_xrot))
elif axis == 'z':
field_rot = np.zeros((3, self.dim[0], self.dim[2], self.dim[1]))
for i in range(self.dim[0]):
mag_x, mag_y, mag_z = self.field[:, i, :, :]
mag_xrot, mag_yrot, mag_zrot = np.rot90(mag_x), np.rot90(mag_y), np.rot90(mag_z)
field_rot[:, i, :, :] = np.array((mag_yrot, -mag_xrot, mag_zrot))
else:
raise ValueError("Wrong input! 'x', 'y', 'z' allowed!")
return VectorData(self.a, field_rot)
def get_slice(self, ax_slice=None, proj_axis='z'):
"""Extract a slice from the :class:`~.VectorData` object.
Parameters
----------
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which the slice is taken. The default is 'z'.
ax_slice : None or int, optional
The slice-index of the axis specified in `proj_axis`. Defaults to the center slice.
Returns
-------
u_mag, v_mag : :class:`~numpy.ndarray` (N=2)
The extracted vector field components in plane perpendicular to the `proj_axis`.
"""
self._log.debug('Calling get_slice')
# Find slice:
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if ax_slice is None:
ax_slice = self.dim[{'z': 0, 'y': 1, 'x': 2}[proj_axis]] // 2
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
self._log.debug('proj_axis == z')
u_mag = np.copy(self.field[0][ax_slice, ...]) # x-component
v_mag = np.copy(self.field[1][ax_slice, ...]) # y-component
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
self._log.debug('proj_axis == y')
u_mag = np.copy(self.field[0][:, ax_slice, :]) # x-component
v_mag = np.copy(self.field[2][:, ax_slice, :]) # z-component
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
self._log.debug('proj_axis == x')
u_mag = np.swapaxes(np.copy(self.field[2][..., ax_slice]), 0, 1) # z-component
v_mag = np.swapaxes(np.copy(self.field[1][..., ax_slice]), 0, 1) # y-component
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(proj_axis))
return u_mag, v_mag
def to_signal(self):
"""Convert :class:`~.VectorData` data into a HyperSpy signal.
Returns
-------
signal: :class:`~hyperspy.signals.Signal`
Representation of the :class:`~.VectorData` object as a HyperSpy Signal.
Notes
-----
This method recquires the hyperspy package!
"""
self._log.debug('Calling to_signal')
signal = super().to_signal()
# Set component axis:
signal.axes_manager[3].name = 'x/y/z-component'
signal.axes_manager[3].units = ''
# Set metadata:
signal.metadata.Signal.title = 'VectorData'
# Return signal:
return signal
def save(self, filename, **kwargs):
"""Saves the VectorData in the specified format.
The function gets the format from the extension:
- hdf5 for HDF5.
- EMD Electron Microscopy Dataset format (also HDF5).
- llg format.
- ovf format.
- npy or npz for numpy formats.
If no extension is provided, 'hdf5' is used. Most formats are
saved with the HyperSpy package (internally the fielddata is first
converted to a HyperSpy Signal.
Each format accepts a different set of parameters. For details
see the specific format documentation.
Parameters
----------
filename : str, optional
Name of the file which the VectorData is saved into. The extension
determines the saving procedure.
"""
from .file_io.io_vectordata import save_vectordata
save_vectordata(self, filename, **kwargs)
def plot_field(self, title='Vector Field', axis=None, proj_axis='z', figsize=(8.5, 7),
ax_slice=None, show_mask=True, bgcolor='white', hue_mode='triadic'):
"""Plot a slice of the vector field as a quiver plot.
Parameters
----------
title : string, optional
The title for the plot.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which a slice is plotted. The default is 'z'.
figsize : tuple of floats (N=2)
Size of the plot figure.
ax_slice : int, optional
The slice-index of the axis specified in `proj_axis`. Is set to the center of
`proj_axis` if not specified.
show_mask: boolean
Default is True. Shows the outlines of the mask slice if available.
bgcolor: {'black', 'white'}, optional
Determines the background color of the plot.
hue_mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme. Use either a triadic or tetradic
scheme (see the according colormaps for more information).
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
"""
self._log.debug('Calling plot_field')
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if ax_slice is None:
ax_slice = self.dim[{'z': 0, 'y': 1, 'x': 2}[proj_axis]] // 2
u_mag, v_mag = self.get_slice(ax_slice, proj_axis)
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
u_label = 'x [px]'
v_label = 'y [px]'
submask = self.get_mask()[ax_slice, ...]
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
u_label = 'x [px]'
v_label = 'z [px]'
submask = self.get_mask()[:, ax_slice, :]
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
u_label = 'z [px]'
v_label = 'y [px]'
submask = self.get_mask()[..., ax_slice]
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(proj_axis))
# If no axis is specified, a new figure is created:
if axis is None:
self._log.debug('axis is None')
fig = plt.figure(figsize=figsize)
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the field:
dim_uv = u_mag.shape
hue = np.arctan2(v_mag, u_mag) / (2 * np.pi) # Hue according to angle!
hue[hue < 0] += 1 # Shift negative values!
luminance = 0.5 * submask # Luminance according to mask!
if bgcolor == 'white': # Invert luminance:
luminance = 1 - luminance
saturation = np.hypot(u_mag, v_mag) # Saturation according to amplitude!
saturation /= saturation.max()
rgb = colors.rgb_from_hls(hue, luminance, saturation, mode=hue_mode)
axis.imshow(Image.fromarray(rgb), origin='lower', interpolation='none',
extent=(0, dim_uv[1], 0, dim_uv[0]))
# Change background color:
axis.set_axis_bgcolor(bgcolor)
# Show mask:
if show_mask and not np.all(submask): # Plot mask if desired and not trivial!
vv, uu = np.indices(dim_uv) + 0.5 # shift to center of pixel
mask_color = 'white' if bgcolor == 'black' else 'black'
axis.contour(uu, vv, submask, levels=[0.5], colors=mask_color,
linestyles='dotted', linewidths=2)
# Further plot formatting:
axis.set_xlim(0, dim_uv[1])
axis.set_ylim(0, dim_uv[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
if dim_uv[0] >= dim_uv[1]:
u_bin, v_bin = np.max((2, np.floor(9 * dim_uv[1] / dim_uv[0]))), 9
else:
u_bin, v_bin = 9, np.max((2, np.floor(9 * dim_uv[0] / dim_uv[1])))
axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
# Return plotting axis:
return axis
def plot_streamline(self, title='Vector Field', axis=None, proj_axis='z', figsize=(8.5, 7),
coloring='angle', ax_slice=None, density=2, linewidth=2,
show_mask=True, bgcolor='white', hue_mode='triadic'):
"""Plot a slice of the vector field as a quiver plot.
Parameters
----------
title : string, optional
The title for the plot.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which a slice is plotted. The default is 'z'.
figsize : tuple of floats (N=2)
Size of the plot figure.
coloring : {'angle', 'amplitude', 'uniform'}
Color coding mode of the arrows. Use 'full' (default), 'angle', 'amplitude' or
'uniform'.
ax_slice : int, optional
The slice-index of the axis specified in `proj_axis`. Is set to the center of
`proj_axis` if not specified.
density : float or 2-tuple, optional
Controls the closeness of streamlines. When density = 1, the domain is divided into a
30x30 grid—density linearly scales this grid. Each cebll in the grid can have, at most,
one traversing streamline. For different densities in each direction, use
[density_x, density_y].
linewidth : numeric or 2d array, optional
Vary linewidth when given a 2d array with the same shape as velocities.
show_mask: boolean
Default is True. Shows the outlines of the mask slice if available.
bgcolor: {'black', 'white'}, optional
Determines the background color of the plot.
hue_mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme. Use either a triadic or tetradic
scheme (see the according colormaps for more information).
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
"""
self._log.debug('Calling plot_quiver')
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if ax_slice is None:
ax_slice = self.dim[{'z': 0, 'y': 1, 'x': 2}[proj_axis]] // 2
u_mag, v_mag = self.get_slice(ax_slice, proj_axis)
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
u_label = 'x [px]'
v_label = 'y [px]'
submask = self.get_mask()[ax_slice, ...]
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
u_label = 'x [px]'
v_label = 'z [px]'
submask = self.get_mask()[:, ax_slice, :]
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
u_label = 'z [px]'
v_label = 'y [px]'
submask = self.get_mask()[..., ax_slice]
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(proj_axis))
# Prepare quiver (select only used arrows if ar_dens is specified):
dim_uv = u_mag.shape
uu = np.arange(dim_uv[1]) + 0.5 # shift to center of pixel
vv = np.arange(dim_uv[0]) + 0.5 # shift to center of pixel
u_mag, v_mag = self.get_slice(ax_slice, proj_axis)
# v_mag = np.ma.array(v_mag, mask=submask)
amplitudes = np.hypot(u_mag, v_mag)
# Calculate the arrow colors:
if coloring == 'angle':
self._log.debug('Encoding angles')
color = np.arctan2(v_mag, u_mag) / (2 * np.pi)
color[color < 0] += 1
if hue_mode == 'triadic':
cmap = colors.hls_triadic_cmap
elif hue_mode == 'tetradic':
cmap = colors.hls_tetradic_cmap
else:
raise ValueError('Hue mode {} not understood!'.format(hue_mode))
elif coloring == 'amplitude':
self._log.debug('Encoding amplitude')
color = amplitudes / amplitudes.max()
cmap = 'jet'
elif coloring == 'uniform':
self._log.debug('No color encoding')
color = np.zeros_like(u_mag) # use black arrows!
cmap = 'gray' if bgcolor == 'white' else 'Greys'
else:
raise AttributeError("Invalid coloring mode! Use 'angles', 'amplitude' or 'uniform'!")
# If no axis is specified, a new figure is created:
if axis is None:
self._log.debug('axis is None')
fig = plt.figure(figsize=figsize)
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the streamlines:
im = plt.streamplot(uu, vv, u_mag, v_mag, density=density, linewidth=linewidth,
color=color, cmap=cmap)
if coloring == 'amplitude':
fig = plt.gcf()
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.82, 0.15, 0.02, 0.7])
cbar = fig.colorbar(im.lines, cax=cbar_ax)
cbar.ax.tick_params(labelsize=14)
cbar_title = u'amplitude'
cbar.set_label(cbar_title, fontsize=15)
# Change background color:
axis.set_axis_bgcolor(bgcolor)
# Show mask:
if show_mask and not np.all(submask): # Plot mask if desired and not trivial!
vv, uu = np.indices(dim_uv) + 0.5 # shift to center of pixel
mask_color = 'white' if bgcolor == 'black' else 'black'
axis.contour(uu, vv, submask, levels=[0.5], colors=mask_color,
linestyles='dotted', linewidths=2)
# Further plot formatting:
axis.set_xlim(0, dim_uv[1])
axis.set_ylim(0, dim_uv[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
if dim_uv[0] >= dim_uv[1]:
u_bin, v_bin = np.max((2, np.floor(9 * dim_uv[1] / dim_uv[0]))), 9
else:
u_bin, v_bin = 9, np.max((2, np.floor(9 * dim_uv[0] / dim_uv[1])))
axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
# Return plotting axis:
return axis
def plot_quiver(self, title='Vector Field', axis=None, proj_axis='z', figsize=(8.5, 7),
coloring='angle', ar_dens=1, ax_slice=None, log=False, scaled=True,
scale=1., show_mask=True, bgcolor='white', hue_mode='triadic'):
"""Plot a slice of the vector field as a quiver plot.
Parameters
----------
title : string, optional
The title for the plot.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which a slice is plotted. The default is 'z'.
figsize : tuple of floats (N=2)
Size of the plot figure.
coloring : {'angle', 'amplitude', 'uniform'}
Color coding mode of the arrows. Use 'full' (default), 'angle', 'amplitude' or
'uniform'.
ar_dens: int, optional
Number defining the arrow density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
ax_slice : int, optional
The slice-index of the axis specified in `proj_axis`. Is set to the center of
`proj_axis` if not specified.
log : boolean, optional
The loratihm of the arrow length is plotted instead. This is helpful if only the
direction of the arrows is important and the amplitude varies a lot. Default is False.
scaled : boolean, optional
Normalizes the plotted arrows in respect to the highest one. Default is True.
scale: float, optional
Additional multiplicative factor scaling the arrow length. Default is 1
(no further scaling).
show_mask: boolean
Default is True. Shows the outlines of the mask slice if available.
bgcolor: {'black', 'white'}, optional
Determines the background color of the plot.
hue_mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme. Use either a triadic or tetradic
scheme (see the according colormaps for more information).
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
"""
self._log.debug('Calling plot_quiver')
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if ax_slice is None:
ax_slice = self.dim[{'z': 0, 'y': 1, 'x': 2}[proj_axis]] // 2
u_mag, v_mag = self.get_slice(ax_slice, proj_axis)
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
u_label = 'x [px]'
v_label = 'y [px]'
submask = self.get_mask()[ax_slice, ...]
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
u_label = 'x [px]'
v_label = 'z [px]'
submask = self.get_mask()[:, ax_slice, :]
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
u_label = 'z [px]'
v_label = 'y [px]'
submask = self.get_mask()[..., ax_slice]
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(proj_axis))
# Prepare quiver (select only used arrows if ar_dens is specified):
dim_uv = u_mag.shape
vv, uu = np.indices(dim_uv) + 0.5 # shift to center of pixel
uu = uu[::ar_dens, ::ar_dens]
vv = vv[::ar_dens, ::ar_dens]
u_mag = u_mag[::ar_dens, ::ar_dens]
v_mag = v_mag[::ar_dens, ::ar_dens]
amplitudes = np.hypot(u_mag, v_mag)
angles = np.angle(u_mag + 1j * v_mag, deg=True).tolist()
# Calculate the arrow colors:
if coloring == 'angle':
self._log.debug('Encoding angles')
hue = np.arctan2(v_mag, u_mag) / (2 * np.pi)
hue[hue < 0] += 1
if hue_mode == 'triadic':
cmap = colors.hls_triadic_cmap
elif hue_mode == 'tetradic':
cmap = colors.hls_tetradic_cmap
else:
raise ValueError('Hue mode {} not understood!'.format(hue_mode))
elif coloring == 'amplitude':
self._log.debug('Encoding amplitude')
hue = amplitudes / amplitudes.max()
cmap = 'jet'
elif coloring == 'uniform':
self._log.debug('No color encoding')
hue = np.zeros_like(u_mag) # use black arrows!
cmap = 'gray' if bgcolor == 'white' else 'Greys'
else:
raise AttributeError("Invalid coloring mode! Use 'angles', 'amplitude' or 'uniform'!")
# If no axis is specified, a new figure is created:
if axis is None:
self._log.debug('axis is None')
fig = plt.figure(figsize=figsize)
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Take the logarithm of the arrows to clearly show directions (if specified):
if log and np.any(amplitudes): # If the slice is empty, skip!
cutoff = 10
amp = np.round(amplitudes, decimals=cutoff)
min_value = amp[np.nonzero(amp)].min()
u_mag = np.round(u_mag, decimals=cutoff) / min_value
u_mag = np.log10(np.abs(u_mag) + 1) * np.sign(u_mag)
v_mag = np.round(v_mag, decimals=cutoff) / min_value
v_mag = np.log10(np.abs(v_mag) + 1) * np.sign(v_mag)
amplitudes = np.hypot(u_mag, v_mag) # Recalculate (used if scaled)!
# Scale the amplitude of the arrows to the highest one (if specified):
if scaled:
u_mag /= amplitudes.max() + 1E-30
v_mag /= amplitudes.max() + 1E-30
im = axis.quiver(uu, vv, u_mag, v_mag, hue, cmap=cmap, clim=(0, 1), angles=angles,
pivot='middle', units='xy', scale_units='xy', scale=scale / ar_dens,
minlength=0.05, width=1*ar_dens, headlength=2, headaxislength=2,
headwidth=2, minshaft=2)
if coloring == 'amplitude':
fig = plt.gcf()
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.82, 0.15, 0.02, 0.7])
cbar = fig.colorbar(im, cax=cbar_ax)
cbar.ax.tick_params(labelsize=14)
cbar_title = u'amplitude'
cbar.set_label(cbar_title, fontsize=15)
# Change background color:
axis.set_axis_bgcolor(bgcolor)
# Show mask:
if show_mask and not np.all(submask): # Plot mask if desired and not trivial!
vv, uu = np.indices(dim_uv) + 0.5 # shift to center of pixel
mask_color = 'white' if bgcolor == 'black' else 'black'
axis.contour(uu, vv, submask, levels=[0.5], colors=mask_color,
linestyles='dotted', linewidths=2)
# Further plot formatting:
axis.set_xlim(0, dim_uv[1])
axis.set_ylim(0, dim_uv[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
if dim_uv[0] >= dim_uv[1]:
u_bin, v_bin = np.max((2, np.floor(9 * dim_uv[1] / dim_uv[0]))), 9
else:
u_bin, v_bin = 9, np.max((2, np.floor(9 * dim_uv[0] / dim_uv[1])))
axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
# Return plotting axis:
return axis
def plot_quiver3d(self, title='Vector Field', limit=None, cmap='jet', mode='2darrow',
coloring='full', ar_dens=1, opacity=1.0, hue_mode='triadic'):
"""Plot the vector field as 3D-vectors in a quiverplot.
Parameters
----------
title : string, optional
The title for the plot.
limit : float, optional
Plotlimit for the vector field arrow length used to scale the colormap.
cmap : string, optional
String describing the colormap which is used (default is 'jet').
ar_dens: int, optional
Number defining the arrow density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
mode: string, optional
Mode, determining the glyphs used in the 3D plot. Default is '2darrow', which
corresponds to 2D arrows. For smaller amounts of arrows, 'arrow' (3D) is prettier.
coloring : {'full', 'angle', 'amplitude'}, optional
Color coding mode of the arrows. Use 'angle' (default) or 'amplitude'.
opacity: float, optional
Defines the opacity of the arrows. Default is 1.0 (completely opaque).
hue_mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme. Use either a triadic or tetradic
scheme (see the according colormaps for more information).
Returns
-------
plot : :class:`mayavi.modules.vectors.Vectors`
The plot object.
"""
self._log.debug('Calling quiver_plot3D')
from mayavi import mlab
a = self.a
dim = self.dim
if limit is None:
limit = np.max(np.nan_to_num(self.field_amp))
ad = ar_dens
# Create points and vector components as lists:
zzz, yyy, xxx = (np.indices(dim) - a / 2).reshape((3,) + dim)
zzz = zzz[::ad, ::ad, ::ad].ravel()
yyy = yyy[::ad, ::ad, ::ad].ravel()
xxx = xxx[::ad, ::ad, ::ad].ravel()
x_mag = self.field[0][::ad, ::ad, ::ad].ravel()
y_mag = self.field[1][::ad, ::ad, ::ad].ravel()
z_mag = self.field[2][::ad, ::ad, ::ad].ravel()
# Plot them as vectors:
mlab.figure(size=(750, 700))
if coloring in ('full', 'angle'): # Encodes the full angle via colorwheel and saturation:
self._log.debug('Encoding full 3D angles')
vecs = mlab.quiver3d(xxx, yyy, zzz, x_mag, y_mag, z_mag, mode=mode, opacity=opacity,
scalars=np.arange(len(xxx)))
h, l, s = colors.hls_from_vector(x_mag, y_mag, z_mag)
if coloring == 'angle': # Encode just the angle and not the amplitude via saturation:
s = np.ones_like(s)
rgbs = colors.rgb_from_hls(h, l, s, mode=hue_mode)
rgbas = np.hstack((rgbs, 255 * np.ones((len(xxx), 1)))).astype(np.uint8)
vecs.glyph.color_mode = 'color_by_scalar'
vecs.module_manager.scalar_lut_manager.lut.table = rgbas
mlab.draw()
elif coloring == 'amplitude': # Encodes the amplitude of the arrows with the jet colormap:
self._log.debug('Encoding amplitude')
vecs = mlab.quiver3d(xxx, yyy, zzz, x_mag, y_mag, z_mag,
mode=mode, colormap=cmap, opacity=opacity)
mlab.colorbar(label_fmt='%.2f')
mlab.colorbar(orientation='vertical')
else:
raise AttributeError('Coloring mode not supported!')
vecs.glyph.glyph_source.glyph_position = 'center'
vecs.module_manager.vector_lut_manager.data_range = np.array([0, limit])
mlab.outline(vecs)
mlab.axes(vecs)
mlab.title(title, height=0.95, size=0.35)
mlab.orientation_axes()
return vecs
class ScalarData(FieldData):
"""Class for storing scalar field data.
Represents 3-dimensional scalar field distributions which is stored as a 3-dimensional
numpy array in `field`, but which can also be accessed as a vector via `field_vec`.
:class:`~.ScalarData` objects support negation, arithmetic operators (``+``, ``-``, ``*``)
and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers and other
:class:`~.ScalarData` objects, if their dimensions and grid spacings match. It is possible
to load data from HDF5 or LLG (.txt) files or to save the data in these formats.
Plotting methods are also provided.
Attributes
----------
a: float
The grid spacing in nm.
field: :class:`~numpy.ndarray` (N=4)
The scalar field.
"""
_log = logging.getLogger(__name__ + '.ScalarData')
def scale_down(self, n=1):
"""Scale down the field distribution by averaging over two pixels along each axis.
Parameters
----------
n : int, optional
Number of times the field distribution is scaled down. The default is 1.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
Only possible, if each axis length is a power of 2!
"""
self._log.debug('Calling scale_down')
assert n > 0 and isinstance(n, int), 'n must be a positive integer!'
self.a *= 2 ** n
for t in range(n):
# Pad if necessary:
pz, py, px = self.dim[0] % 2, self.dim[1] % 2, self.dim[2] % 2
if pz != 0 or py != 0 or px != 0:
self.field = np.pad(self.field, ((0, pz), (0, py), (0, px)), mode='constant')
# Create coarser grid for the field:
shape_4d = (self.dim[0] / 2, 2, self.dim[1] / 2, 2, self.dim[2] / 2, 2)
self.field = self.field.reshape(shape_4d).mean(axis=(5, 3, 1))
def scale_up(self, n=1, order=0):
"""Scale up the field distribution using spline interpolation of the requested order.
Parameters
----------
n : int, optional
Power of 2 with which the grid is scaled. Default is 1, which means every axis is
increased by a factor of ``2**1 = 2``.
order : int, optional
The order of the spline interpolation, which has to be in the range between 0 and 5
and defaults to 0.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
"""
self._log.debug('Calling scale_up')
assert n > 0 and isinstance(n, int), 'n must be a positive integer!'
assert 5 > order >= 0 and isinstance(order, int), \
'order must be a positive integer between 0 and 5!'
self.a /= 2 ** n
self.field = zoom(self.field, zoom=2 ** n, order=order)
def get_vector(self, mask):
"""Returns the field as a vector, specified by a mask.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Masks the pixels from which the components should be taken.
Returns
-------
vector : :class:`~numpy.ndarray` (N=1)
The vector containing the field of the specified pixels.
"""
self._log.debug('Calling get_vector')
if mask is not None:
return np.reshape(self.field[mask], -1)
else:
return self.field_vec
def set_vector(self, vector, mask=None):
"""Set the field components of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean), optional
Masks the pixels from which the components should be taken.
vector : :class:`~numpy.ndarray` (N=1)
The vector containing the field of the specified pixels.
Returns
-------
None
"""
self._log.debug('Calling set_vector')
vector = np.asarray(vector, dtype=fft.FLOAT)
if mask is not None:
self.field[mask] = vector
else:
self.field_vec = vector
def to_signal(self):
"""Convert :class:`~.ScalarData` data into a HyperSpy signal.
Returns
-------
signal: :class:`~hyperspy.signals.Signal`
Representation of the :class:`~.ScalarData` object as a HyperSpy Signal.
Notes
-----
This method recquires the hyperspy package!
"""
self._log.debug('Calling to_signal')
signal = super().to_signal()
# Set metadata:
signal.metadata.Signal.title = 'ScalarData'
# Return signal:
return signal
def save(self, filename, **kwargs):
"""Saves the ScalarData in the specified format.
The function gets the format from the extension:
- hdf5 for HDF5.
- EMD Electron Microscopy Dataset format (also HDF5).
- npy or npz for numpy formats.
If no extension is provided, 'hdf5' is used. Most formats are
saved with the HyperSpy package (internally the fielddata is first
converted to a HyperSpy Signal.
Each format accepts a different set of parameters. For details
see the specific format documentation.
Parameters
----------
filename : str, optional
Name of the file which the ScalarData is saved into. The extension
determines the saving procedure.
"""
from .file_io.io_scalardata import save_scalardata
save_scalardata(self, filename, **kwargs)
pyramid/file_io/__init__.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing Pyramid IO functionality."""
from .io_phasemap import load_phasemap
from .io_vectordata import load_vectordata
from .io_scalardata import load_scalardata
__all__ = ['load_phasemap', 'load_vectordata', 'load_scalardata']
pyramid/file_io/io_phasemap.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""IO functionality for FieldData objects."""
import logging
import os
import numpy as np
from ..phasemap import PhaseMap
__all__ = ['load_phasemap']
_log = logging.getLogger(__name__)
def load_phasemap(filename, mask=None, confidence=None, a=None, **kwargs):
"""Load supported file into a :class:`~pyramid.phasemap.PhaseMap` instance.
The function loads the file according to the extension:
- hdf5 for HDF5.
- rpl for Ripple (useful to export to Digital Micrograph).
- dm3 and dm4 for Digital Micrograph files.
- unf for SEMPER unf binary format.
- txt format.
- npy or npz for numpy formats.
- Many image formats such as png, tiff, jpeg...
Any extra keyword is passed to the corresponsing reader. For
available options see their individual documentation.
Parameters
----------
filename: str
The filename to be loaded.
mask: str or tuple of str and dict, optional
If this is a filename, a mask will be loaded from an appropriate file. If additional
keywords should be provided, use a tuple of the filename and an according dictionary.
If this is `None`, no mask will be loaded.
default is False.
confidence: str or tuple of str and dict, optional
If this is a filename, a confidence matrix will be loaded from an appropriate file. If
additional keywords should be provided, use a tuple of the filename and an according
dictionary. If this is `None`, no mask will be loaded.
a: float or None, optional
If the grid spacing is not None, it will override a possibly loaded value from the files.
Returns
-------
phasemap : :class:`~.PhaseMap`
A :class:`~.PhaseMap` object containing the loaded data.
"""
_log.debug('Calling load_phasemap')
phasemap = _load(filename, as_phasemap=True, a=a, **kwargs)
if mask is not None:
filemask, kwargs_mask = _parse_add_param(mask)
phasemap.mask = _load(filemask, **kwargs_mask)
if confidence is not None:
fileconf, kwargs_conf = _parse_add_param(confidence)
phasemap.confidence = _load(fileconf, **kwargs_conf)
return phasemap
def _load(filename, as_phasemap=False, a=1., **kwargs):
_log.debug('Calling _load')
extension = os.path.splitext(filename)[1]
# Load from txt-files:
if extension == '.txt':
return _load_from_txt(filename, as_phasemap, a, **kwargs)
# Load from npy-files:
elif extension in ['.npy', '.npz']:
return _load_from_npy(filename, as_phasemap, a, **kwargs)
# Load with HyperSpy:
else:
if extension == '':
filename = '{}.hdf5'.format(filename) # Default!
return _load_from_hs(filename, as_phasemap, a, **kwargs)
def _parse_add_param(param):
if param is None:
return None, {}
elif isinstance(param, str):
return param, {}
elif isinstance(param, (list, tuple)) and len(param) == 2:
return param
else:
raise ValueError('Parameter can be a string or a tuple/list of a string and a dict!')
def _load_from_txt(filename, as_phasemap, a, **kwargs):
def _load_arr(filename, **kwargs):
with open(filename, 'r') as phase_file:
if phase_file.readline().startswith('PYRAMID'): # File has pyramid structure:
return np.loadtxt(filename, delimiter='\t', skiprows=2)
else: # Try default args:
return np.loadtxt(filename, **kwargs)
result = _load_arr(filename, **kwargs)
if as_phasemap:
if a is None:
a = 1. # Default!
with open(filename, 'r') as phase_file:
header = phase_file.readline()
if header.startswith('PYRAMID'): # File has pyramid structure:
a = float(phase_file.readline()[15:-4])
return PhaseMap(a, result)
else:
return result
def _load_from_npy(filename, as_phasemap, a, **kwargs):
result = np.load(filename, **kwargs)
if as_phasemap:
if a is None:
a = 1. # Use default!
return PhaseMap(a, result)
else:
return result
def _load_from_hs(filename, as_phasemap, a, **kwargs):
try:
import hyperspy.api as hs
except ImportError:
_log.error('This method recquires the hyperspy package!')
return
phasemap = PhaseMap.from_signal(hs.load(filename, **kwargs))
if as_phasemap:
if a is not None:
phasemap.a = a
return phasemap
else:
return phasemap.phase
def save_phasemap(phasemap, filename, save_mask, save_conf, pyramid_format, **kwargs):
"""%s"""
_log.debug('Calling save_phasemap')
extension = os.path.splitext(filename)[1]
if extension == '.txt': # Save to txt-files:
_save_to_txt(phasemap, filename, pyramid_format, save_mask, save_conf, **kwargs)
elif extension in ['.npy', '.npz']: # Save to npy-files:
_save_to_npy(phasemap, filename, save_mask, save_conf, **kwargs)
else: # Try HyperSpy:
_save_to_hs(phasemap, filename, save_mask, save_conf, **kwargs)
save_phasemap.__doc__ %= PhaseMap.save.__doc__
def _save_to_txt(phasemap, filename, pyramid_format, save_mask, save_conf, **kwargs):
def _save_arr(filename, array, tag, **kwargs):
if pyramid_format:
with open(filename, 'w') as phase_file:
name = os.path.splitext(os.path.split(filename)[1])[0]
phase_file.write('PYRAMID-{}: {}\n'.format(tag, name))
phase_file.write('grid spacing = {} nm\n'.format(phasemap.a))
save_kwargs = {'fmt': '%7.6e', 'delimiter': '\t'}
else:
save_kwargs = kwargs
with open(filename, 'ba') as phase_file:
np.savetxt(phase_file, array, **save_kwargs)
name, extension = os.path.splitext(filename)
_save_arr('{}{}'.format(name, extension), phasemap.phase, 'PHASEMAP', **kwargs)
if save_mask:
_save_arr('{}_mask{}'.format(name, extension), phasemap.mask, 'MASK', **kwargs)
if save_conf:
_save_arr('{}_conf{}'.format(name, extension), phasemap.confidence, 'CONFIDENCE', **kwargs)
def _save_to_npy(phasemap, filename, save_mask, save_conf, **kwargs):
name, extension = os.path.splitext(filename)
np.save('{}{}'.format(name, extension), phasemap.phase, **kwargs)
if save_mask:
np.save('{}_mask{}'.format(name, extension), phasemap.mask, **kwargs)
if save_conf:
np.save('{}_conf{}'.format(name, extension), phasemap.confidence, **kwargs)
def _save_to_hs(phasemap, filename, save_mask, save_conf, **kwargs):
name, extension = os.path.splitext(filename)
phasemap.to_signal().save(filename, **kwargs)
if extension not in ['.hdf5', '.HDF5', '']: # mask and confidence are saved separately:
import hyperspy.api as hs
if save_mask:
mask_name = '{}_mask{}'.format(name, extension)
hs.signals.Signal2D(phasemap.mask, **kwargs).save(mask_name)
if save_conf:
conf_name = '{}_conf{}'.format(name, extension)
hs.signals.Signal2D(phasemap.confidence, **kwargs).save(conf_name)
pyramid/file_io/io_scalardata.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""IO functionality for ScalarData objects."""
import logging
import os
import numpy as np
from ..fielddata import ScalarData
__all__ = ['load_scalardata']
_log = logging.getLogger(__name__)
def load_scalardata(filename, a=None, **kwargs):
"""Load supported file into a :class:`~pyramid.fielddata.ScalarData` instance.
The function loads the file according to the extension:
- hdf5 for HDF5.
- EMD Electron Microscopy Dataset format (also HDF5).
- npy or npz for numpy formats.
Any extra keyword is passed to the corresponsing reader. For
available options see their individual documentation.
Parameters
----------
filename: str
The filename to be loaded.
a: float or None, optional
If the grid spacing is not None, it will override a possibly loaded value from the files.
Returns
-------
scalardata : :class:`~.ScalarData`
A :class:`~.ScalarData` object containing the loaded data.
"""
_log.debug('Calling load_scalardata')
extension = os.path.splitext(filename)[1]
# Load from npy-files:
if extension in ['.npy', '.npz']:
return _load_from_npy(filename, a, **kwargs)
# Load with HyperSpy:
else:
if extension == '':
filename = '{}.hdf5'.format(filename) # Default!
return _load_from_hs(filename, a, **kwargs)
def _load_from_npy(filename, a, **kwargs):
_log.debug('Calling load_from_npy')
if a is None:
a = 1. # Use default!
return ScalarData(a, np.load(filename, **kwargs))
def _load_from_hs(filename, a, **kwargs):
_log.debug('Calling load_from_hs')
try:
import hyperspy.api as hs
except ImportError:
_log.error('This method recquires the hyperspy package!')
return
scalardata = ScalarData.from_signal(hs.load(filename, **kwargs))
if a is not None:
scalardata.a = a
return scalardata
def save_scalardata(scalardata, filename, **kwargs):
"""%s"""
_log.debug('Calling save_scalardata')
extension = os.path.splitext(filename)[1]
if extension in ['.npy', '.npz']: # Save to npy-files:
_save_to_npy(scalardata, filename, **kwargs)
else: # Try HyperSpy:
_save_to_hs(scalardata, filename, **kwargs)
save_scalardata.__doc__ %= ScalarData.save.__doc__
def _save_to_npy(scalardata, filename, **kwargs):
_log.debug('Calling save_to_npy')
np.save(filename, scalardata.field, **kwargs)
def _save_to_hs(scalardata, filename, **kwargs):
_log.debug('Calling save_to_hyperspy')
scalardata.to_signal().save(filename, **kwargs)
pyramid/file_io/io_vectordata.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""IO functionality for VectorData objects."""
import logging
import os
import numpy as np
from ..fielddata import VectorData
__all__ = ['load_vectordata']
_log = logging.getLogger(__name__)
def load_vectordata(filename, a=None, **kwargs):
"""Load supported file into a :class:`~pyramid.fielddata.VectorData` instance.
The function loads the file according to the extension:
- hdf5 for HDF5.
- EMD Electron Microscopy Dataset format (also HDF5).
- llg format.
- ovf format.
- npy or npz for numpy formats.
Any extra keyword is passed to the corresponsing reader. For
available options see their individual documentation.
Parameters
----------
filename: str
The filename to be loaded.
a: float or None, optional
If the grid spacing is not None, it will override a possibly loaded value from the files.
Returns
-------
vectordata : :class:`~.VectorData`
A :class:`~.VectorData` object containing the loaded data.
"""
_log.debug('Calling load_vectordata')
extension = os.path.splitext(filename)[1]
# Load from llg-files:
if extension in ['.llg', '.txt']:
return _load_from_llg(filename, a)
# Load from ovf-files:
if extension in ['.ovf', '.omf', '.ohf', 'obf']:
return _load_from_ovf(filename, a)
# Load from npy-files:
elif extension in ['.npy', '.npz']:
return _load_from_npy(filename, a, **kwargs)
# Load with HyperSpy:
else:
if extension == '':
filename = '{}.hdf5'.format(filename) # Default!
return _load_from_hs(filename, a, **kwargs)
def _load_from_llg(filename, a):
_log.debug('Calling load_from_llg')
SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
data = np.genfromtxt(filename, skip_header=2)
dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data[:, 0])))
if a is None:
a = (data[1, 0] - data[0, 0]) / SCALE
field = data[:, 3:6].T.reshape((3,) + dim)
return VectorData(a, field)
def _load_from_ovf(filename, a):
_log.debug('Calling load_from_ovf')
with open(filename, 'r') as mag_file:
assert mag_file.readline().startswith('# OOMMF') # Make sure file has .ovf-format!
read_header, read_data = False, False
header = {}
x_mag, y_mag, z_mag = [], [], []
for line in mag_file:
# Read in additional info:
if not read_header and not read_data:
if line.startswith('# Segment count'):
assert int(line.split()[3]) == 1, 'Only one vector field can be read at time!'
elif line.startswith('# Begin: Header'):
read_header = True
elif line.startswith('# Begin: Data'):
read_data = True
data_mode = line.split()[3]
assert data_mode in ['text', 'Text'], \
'Data mode {} is currently not supported by this reader!'.format(data_mode)
assert header.get('meshtype') == 'rectangular', \
'Only rectangular grids can be currently read!'
# Read in header:
elif read_header: # Read header:
if line.startswith('# End: Header'): # Header is done:
read_header = False
continue
line_list = line.split()
if '##' in line_list: # Strip trailing comments:
del line_list[line_list.index('##'):]
if len(line_list) <= 1: # Just '#' or empty line:
continue
key, value = line_list[1].strip(':'), ' '.join(line_list[2:])
if key not in header: # Add new key, value pair:
header[key] = value
elif key == 'Desc': # Can go over several lines:
header['Desc'] = ' '.join([header['Desc'], value])
# Read in data:
elif read_data: # Currently in a data block:
if line.startswith('# End: Data'): # Header is done:
read_data = False
else:
x, y, z = [float(i) for i in line.split()]
x_mag.append(x)
y_mag.append(y)
z_mag.append(z)
# Format after reading:
dim = (int(header['znodes']), int(header['ynodes']), int(header['xnodes']))
x_mag = np.asarray(x_mag).reshape(dim)
y_mag = np.asarray(y_mag).reshape(dim)
z_mag = np.asarray(z_mag).reshape(dim)
field = np.asarray((x_mag, y_mag, z_mag)) * float(header.get('valuemultiplier', 1))
if a is None:
assert header.get('xstepsize') == header.get('ystepsize') == header.get('zstepsize'), \
'Grid spacing is not equidistant!'
a = float(header.get('xstepsize', 1.))
meshunit = header.get('meshunit', 'nm')
a *= {'m': 1e9, 'mm': 1e6, 'µm': 1e3, 'nm': 1}[meshunit] # Conversion to nm
return VectorData(a, field)
def _load_from_npy(filename, a, **kwargs):
_log.debug('Calling load_from_npy')
if a is None:
a = 1. # Use default!
return VectorData(a, np.load(filename, **kwargs))
def _load_from_hs(filename, a, **kwargs):
_log.debug('Calling load_from_hs')
try:
import hyperspy.api as hs
except ImportError:
_log.error('This method recquires the hyperspy package!')
return
vectordata = VectorData.from_signal(hs.load(filename, **kwargs))
if a is not None:
vectordata.a = a
return vectordata
def save_vectordata(vectordata, filename, **kwargs):
"""%s"""
_log.debug('Calling save_vectordata')
extension = os.path.splitext(filename)[1]
if extension == '.llg': # Save to llg-files:
_save_to_llg(vectordata, filename)
elif extension == '.ovf': # Save to ovf-files:
_save_to_ovf(vectordata, filename, **kwargs)
elif extension in ['.npy', '.npz']: # Save to npy-files:
_save_to_npy(vectordata, filename, **kwargs)
else: # Try HyperSpy:
_save_to_hs(vectordata, filename, **kwargs)
save_vectordata.__doc__ %= VectorData.save.__doc__
def _save_to_llg(vectordata, filename):
_log.debug('Calling save_to_llg')
dim = vectordata.dim
SCALE = 1.0E-9 / 1.0E-2 # from nm to cm
# Create 3D meshgrid and reshape it and the field into a list where x varies first:
zz, yy, xx = vectordata.a * SCALE * (np.indices(dim) + 0.5).reshape(3, -1)
x_vec, y_vec, z_vec = vectordata.field.reshape(3, -1)
data = np.array([xx, yy, zz, x_vec, y_vec, z_vec]).T
# Save data to file:
with open(filename, 'w') as mag_file:
mag_file.write('LLGFileCreator: %s\n'.format(filename))
mag_file.write(' %d %d %d\n'.format(np.asarray(dim)[::-1]))
mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
for cell in row) for row in data))
def _save_to_ovf(vectordata, filename):
_log.debug('Calling save_to_ovf')
with open(filename, 'w') as mag_file:
mag_file.write('# OOMMF OVF 2.0\n')
mag_file.write('# Segment count: 1\n')
mag_file.write('# Begin: Segment\n')
# Write Header:
mag_file.write('# Begin: Header\n')
name = os.path.split(os.path.split(filename)[1])
mag_file.write('# Title: PYRAMID-VECTORDATA {}\n'.format(name))
mag_file.write('# meshtype: rectangular\n')
mag_file.write('# meshunit: nm\n')
mag_file.write('# valueunit: A/m\n')
mag_file.write('# valuemultiplier: 1.\n')
mag_file.write('# xmin: 0.\n')
mag_file.write('# ymin: 0.\n')
mag_file.write('# zmin: 0.\n')
mag_file.write('# xmax: {}\n'.format(vectordata.a * vectordata.dim[2]))
mag_file.write('# ymax: {}\n'.format(vectordata.a * vectordata.dim[1]))
mag_file.write('# zmax: {}\n'.format(vectordata.a * vectordata.dim[0]))
mag_file.write('# xbase: 0.\n')
mag_file.write('# ybase: 0.\n')
mag_file.write('# zbase: 0.\n')
mag_file.write('# xstepsize: {}\n'.format(vectordata.a))
mag_file.write('# ystepsize: {}\n'.format(vectordata.a))
mag_file.write('# zstepsize: {}\n'.format(vectordata.a))
mag_file.write('# xnodes: {}\n'.format(vectordata.dim[2]))
mag_file.write('# ynodes: {}\n'.format(vectordata.dim[1]))
mag_file.write('# znodes: {}\n'.format(vectordata.dim[0]))
mag_file.write('# End: Header\n')
# Write data:
mag_file.write('# Begin: Data Text\n')
x_mag, y_mag, z_mag = vectordata.field
x_mag = x_mag.ravel()
y_mag = y_mag.ravel()
z_mag = z_mag.ravel()
for i in range(np.prod(vectordata.dim)):
mag_file.write('{:g} {:g} {:g}\n'.format(x_mag[i], y_mag[i], z_mag[i]))
mag_file.write('# End: Data Text\n')
mag_file.write('# End: Segment\n')
def _save_to_npy(vectordata, filename, **kwargs):
_log.debug('Calling save_to_npy')
np.save(filename, vectordata.field, **kwargs)
def _save_to_hs(vectordata, filename, **kwargs):
_log.debug('Calling save_to_hyperspy')
vectordata.to_signal().save(filename, **kwargs)
pyramid/forwardmodel.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.ForwardModel` class which represents a strategy to map a
threedimensional magnetization distribution onto a two-dimensional phase map."""
import logging
import multiprocessing as mp
import sys
import numpy as np
from pyramid.dataset import DataSet
from pyramid.fielddata import VectorData
from pyramid.ramp import Ramp
__all__ = ['ForwardModel', 'DistributedForwardModel']
class ForwardModel(object):
"""Class for mapping 3D magnetic distributions to 2D phase maps.
Represents a strategy for the mapping of a 3D magnetic distribution to two-dimensional
phase maps. A :class:`~.DataSet` object is given which is used as input for the model
(projectors, phasemappers, etc.). A `ramp_order` can be specified to add polynomial ramps
to the constructed phase maps (which can also be reconstructed!). A :class:`~.Ramp` class
object will be constructed accordingly, which also holds all info about the ramps after a
reconstruction.
Attributes
----------
data_set: :class:`~dataset.DataSet`
:class:`~dataset.DataSet` object, which stores all required information calculation.
ramp_order : int or None (default)
Polynomial order of the additional phase ramp which will be added to the phase maps.
All ramp parameters have to be at the end of the input vector and are split automatically.
Default is None (no ramps are added).
y : :class:`~numpy.ndarray` (N=1)
Vector which lists all pixel values of all phase maps one after another.
m: int
Size of the image space. Number of pixels of the 2-dimensional projected grid.
n: int
Size of the input space. Number of voxels of the 3-dimensional grid.
Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `m x m` with m
being the length of the targetvector y (vectorized phase map information).
"""
_log = logging.getLogger(__name__ + '.ForwardModel')
def __init__(self, data_set, ramp_order=None):
self._log.debug('Calling __init__')
self.data_set = data_set
self.ramp_order = ramp_order
# Extract information from data_set:
self.phasemappers = self.data_set.phasemappers
self.y = self.data_set.phase_vec
self.n = self.data_set.n
self.m = self.data_set.m
self.shape = (self.m, self.n)
self.hook_points = self.data_set.hook_points
if self.data_set.Se_inv is None:
self.data_set.set_Se_inv_diag_with_conf()
self.Se_inv = self.data_set.Se_inv
# Create ramp and change n accordingly:
self.ramp = Ramp(self.data_set, self.ramp_order)
self.n += self.ramp.n # ramp.n is 0 if ramp_order is None
# Create empty MagData object:
self.magdata = VectorData(self.data_set.a, np.zeros((3,) + self.data_set.dim))
self._log.debug('Creating ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(data_set=%r)' % (self.__class__, self.data_set)
def __str__(self):
self._log.debug('Calling __str__')
return 'ForwardModel(data_set=%s)' % self.data_set
def __call__(self, x):
# Extract ramp parameters if necessary (x will be shortened!):
x = self.ramp.extract_ramp_params(x)
# Reset magdata and fill with vector:
self.magdata.field[...] = 0
self.magdata.set_vector(x, self.data_set.mask)
# Simulate all phase maps and create result vector:
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
mapper = self.phasemappers[projector.dim_uv]
phasemap = mapper(projector(self.magdata))
phasemap += self.ramp(i) # add ramp!
result[hp[i]:hp[i + 1]] = phasemap.phase_vec
return np.reshape(result, -1)
def jac_dot(self, x, vector):
"""Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Evaluation point of the jacobi-matrix. The Jacobi matrix is constant for a linear
problem, thus `x` can be set to None (it is not used int the computation). It is
implemented for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the 3D magnetization distribution. First the `x`, then the `y` and
lastly the `z` components are listed. Ramp parameters are also added at the end if
necessary.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the input
`vector`.
"""
# Extract ramp parameters if necessary (vector will be shortened!):
vector = self.ramp.extract_ramp_params(vector)
# Reset magdata and fill with vector:
self.magdata.field[...] = 0
self.magdata.set_vector(vector, self.data_set.mask)
# Simulate all phase maps and create result vector:
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
mag_vec = self.magdata.field_vec
mapper = self.phasemappers[projector.dim_uv]
res = mapper.jac_dot(projector.jac_dot(mag_vec))
res += self.ramp.jac_dot(i) # add ramp!
result[hp[i]:hp[i + 1]] = res
return result
def jac_T_dot(self, x, vector):
"""'Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Evaluation point of the jacobi-matrix. The jacobi matrix is constant for a linear
problem, thus `x` can be set to None (it is not used int the computation). Is used
for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of all 2D phase maps one after another in one vector.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitely calculated) with
the input `vector`. If necessary, transposed ramp parameters are concatenated.
"""
proj_T_result = np.zeros(3 * np.prod(self.data_set.dim))
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
sub_vec = vector[hp[i]:hp[i + 1]]
mapper = self.phasemappers[projector.dim_uv]
proj_T_result += projector.jac_T_dot(mapper.jac_T_dot(sub_vec))
self.magdata.field_vec = proj_T_result
result = self.magdata.get_vector(self.data_set.mask)
ramp_params = self.ramp.jac_T_dot(vector) # calculate ramp_params separately!
return np.concatenate((result, ramp_params))
def finalize(self):
"""'Finalize the processes and let them join the master process (NOT USED HERE!).
Returns
-------
None
"""
pass
class DistributedForwardModel(ForwardModel):
"""Multiprocessing class for mapping 3D magnetic distributions to 2D phase maps.
Subclass of the :class:`~.ForwardModel` class which implements multiprocessing strategies
to speed up the calculations. The interface is the same, internally, the processes and one
ForwardModel operating on a subset of the DataSet per process are created during construction.
Ramps are calculated in the main thread. The :func:`~.finalize` method can be used to force
the processes to join if the class is no longer used.
Attributes
----------
data_set: :class:`~dataset.DataSet`
:class:`~dataset.DataSet` object, which stores all required information calculation.
ramp_order : int or None (default)
Polynomial order of the additional phase ramp which will be added to the phase maps.
All ramp parameters have to be at the end of the input vector and are split automatically.
Default is None (no ramps are added).
nprocs: int
Number of processes which should be created. Default is 1 (not recommended).
"""
def __init__(self, data_set, ramp_order=None, nprocs=1):
# Evoke super constructor to set up the normal ForwardModel:
super().__init__(data_set, ramp_order)
# Initialize multirocessing specific stuff:
self.nprocs = nprocs
img_per_proc = np.ceil(self.data_set.count / self.nprocs).astype(np.int)
hp = self.data_set.hook_points
self.proc_hook_points = [0]
self.pipes = []
self.processes = []
for proc_id in range(self.nprocs):
# Create SubDataSets:
sub_data = DataSet(self.data_set.a, self.data_set.dim, self.data_set.b_0,
self.data_set.mask, self.data_set.Se_inv)
# Distribute data to SubDataSets:
start = proc_id * img_per_proc
stop = np.min(((proc_id + 1) * img_per_proc, self.data_set.count))
self.proc_hook_points.append(hp[stop])
sub_data.phasemaps = self.data_set.phasemaps[start:stop]
sub_data.projectors = self.data_set.projectors[start:stop]
# Create SubForwardModel:
sub_fwd_model = ForwardModel(sub_data, ramp_order=None) # ramps handled in master!
# Create communication pipe:
self.pipes.append(mp.Pipe())
# Create process:
p = mp.Process(name='Worker {}'.format(proc_id), target=_worker,
args=(sub_fwd_model, self.pipes[proc_id][1]))
self.processes.append(p)
# Start process:
p.start()
self._log.debug('Creating ' + str(self))
def __call__(self, x):
# Extract ramp parameters if necessary (x will be shortened!):
x = self.ramp.extract_ramp_params(x)
# Distribute input to processes and start working:
for proc_id in range(self.nprocs):
self.pipes[proc_id][0].send(('__call__', (x,)))
# Initialize result vector and shorten hook point names:
result = np.zeros(self.m)
hp = self.hook_points
php = self.proc_hook_points
# Calculate ramps (if necessary):
if self.ramp_order is not None:
for i in range(self.data_set.count):
result[hp[i]:hp[i + 1]] += self.ramp(i).phase.ravel()
# Get process results from the pipes:
for proc_id in range(self.nprocs):
result[php[proc_id]:php[proc_id + 1]] += self.pipes[proc_id][0].recv()
# Return result:
return result
def jac_dot(self, x, vector):
"""Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Evaluation point of the jacobi-matrix. The Jacobi matrix is constant for a linear
problem, thus `x` can be set to None (it is not used int the computation). It is
implemented for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the 3D magnetization distribution. First the `x`, then the `y` and
lastly the `z` components are listed. Ramp parameters are also added at the end if
necessary.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the input
`vector`.
"""
# Extract ramp parameters if necessary (x will be shortened!):
vector = self.ramp.extract_ramp_params(vector)
# Distribute input to processes and start working:
for proc_id in range(self.nprocs):
self.pipes[proc_id][0].send(('jac_dot', (None, vector)))
# Initialize result vector and shorten hook point names:
result = np.zeros(self.m)
hp = self.hook_points
php = self.proc_hook_points
# Calculate ramps (if necessary):
if self.ramp_order is not None:
for i in range(self.data_set.count):
result[hp[i]:hp[i + 1]] += self.ramp.jac_dot(i)
# Get process results from the pipes:
for proc_id in range(self.nprocs):
result[php[proc_id]:php[proc_id + 1]] += self.pipes[proc_id][0].recv()
# Return result:
return result
def jac_T_dot(self, x, vector):
"""'Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Evaluation point of the jacobi-matrix. The jacobi matrix is constant for a linear
problem, thus `x` can be set to None (it is not used int the computation). Is used
for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of all 2D phase maps one after another in one vector.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitely calculated) with
the input `vector`. If necessary, transposed ramp parameters are concatenated.
"""
php = self.proc_hook_points
# Distribute input to processes and start working:
for proc_id in range(self.nprocs):
sub_vec = vector[php[proc_id]:php[proc_id + 1]]
self.pipes[proc_id][0].send(('jac_T_dot', (None, sub_vec)))
# Calculate ramps:
ramp_params = self.ramp.jac_T_dot(vector) # calculate ramp_params separately!
# Initialize result vector:
result = np.zeros(3 * self.data_set.mask.sum())
# Get process results from the pipes:
for proc_id in range(self.nprocs):
result += self.pipes[proc_id][0].recv()
# Return result:
return np.concatenate((result, ramp_params))
def finalize(self):
"""'Finalize the processes and let them join the master process.
Returns
-------
None
"""
# Finalize processes:
for proc_id in range(self.nprocs):
self.pipes[proc_id][0].send('STOP')
# Exit the completed processes:
for p in self.processes:
p.join()
def _worker(fwd_model, pipe):
# Has to be directly accessible in the module as a function, NOT a method of a class instance!
print('... {} starting!'.format(mp.current_process().name))
sys.stdout.flush()
for method, arguments in iter(pipe.recv, 'STOP'):
# '... {} processes method {}'.format(mp.current_process().name, method)
sys.stdout.flush()
result = getattr(fwd_model, method)(*arguments)
pipe.send(result)
print('... ', mp.current_process().name, 'exiting!')
sys.stdout.flush()
pyramid/kernel.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.Kernel` class, representing the phase contribution of one
single magnetized pixel."""
import logging
import numpy as np
from pyramid import fft
__all__ = ['Kernel', 'PHI_0']
PHI_0 = 2067.83 # magnetic flux in T*nm²
class Kernel(object):
"""Class for calculating kernel matrices for the phase calculation.
Represents the phase of a single magnetized pixel for two orthogonal directions (`u` and `v`),
which can be accessed via the corresponding attributes. The default elementary geometry is
`disc`, but can also be specified as the phase of a `slab` representation of a single
magnetized pixel. During the construction, a few attributes are calculated that are used in
the convolution during phase calculation in the different :class:`~Phasemapper` classes.
An instance of the :class:`~.Kernel` class can be called as a function with a `vector`,
which represents the projected magnetization onto a 2-dimensional grid.
Attributes
----------
a : float
The grid spacing in nm.
dim_uv : tuple of int (N=2), optional
Dimensions of the 2-dimensional projected magnetization grid from which the phase should
be calculated.
dim_kern : tuple of int (N=2)
Dimensions of the kernel, which is ``2N-1`` for both axes compared to `dim_uv`.
dim_pad : tuple of int (N=2)
Dimensions of the padded FOV, which is ``2N`` (if FFTW is used) or the next highest power
of 2 (for numpy-FFT).
dim_fft : tuple of int (N=2)
Dimensions of the grid, which is used for the FFT, taking into account that a RFFT should
be used (one axis is halved in comparison to `dim_pad`).
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
u : :class:`~numpy.ndarray` (N=3)
The phase contribution of one pixel magnetized in u-direction.
v : :class:`~numpy.ndarray` (N=3)
The phase contribution of one pixel magnetized in v-direction.
u_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the phase contribution of one pixel magnetized in u-direction.
v_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the phase contribution of one pixel magnetized in v-direction.
slice_phase : tuple (N=2) of :class:`slice`
A tuple of :class:`slice` objects to extract the original FOV from the increased one with
size `dim_pad` for the elementary kernel phase. The kernel is shifted, thus the center is
not at (0, 0), which also shifts the slicing compared to `slice_mag`.
slice_mag : tuple (N=2) of :class:`slice`
A tuple of :class:`slice` objects to extract the original FOV from the increased one with
size `dim_pad` for the projected magnetization distribution.
prw_vec: tuple of 2 int, optional
A two-component vector describing the displacement of the reference wave to include
perturbation of this reference by the object itself (via fringing fields).
"""
_log = logging.getLogger(__name__ + '.Kernel')
def __init__(self, a, dim_uv, b_0=1., prw_vec=None, geometry='disc'):
self._log.debug('Calling __init__')
# Set basic properties:
self.dim_uv = dim_uv # Dimensions of the FOV
self.dim_kern = tuple(2 * np.array(dim_uv) - 1) # Dimensions of the kernel
self.a = a
self.geometry = geometry
# Set up FFT:
if fft.BACKEND == 'pyfftw':
self.dim_pad = tuple(2 * np.array(dim_uv)) # is at least even (not nec. power of 2)
elif fft.BACKEND == 'numpy':
self.dim_pad = tuple(2 ** np.ceil(np.log2(2 * np.array(dim_uv))).astype(int)) # pow(2)
self.dim_fft = (self.dim_pad[0], self.dim_pad[1] // 2 + 1) # last axis is real
self.slice_phase = (slice(dim_uv[0] - 1, self.dim_kern[0]), # Shift because kernel center
slice(dim_uv[1] - 1, self.dim_kern[1])) # is not at (0, 0)!
self.slice_mag = (slice(0, dim_uv[0]), # Magnetization is padded on the far end!
slice(0, dim_uv[1])) # (Phase cutout is shifted as listed above)
# Calculate kernel (single pixel phase):
# [M_0] = A/m --> This is the magnetization, not the magnetic moment (A/m * m³ = Am²)!
# [PHI_0 / µ_0] = Tm² / Tm/A = Am
# [b_0] = [M_0] * [µ_0] = A/m * N/A² = N/Am = T
# [coeff] = [b_0 * a² / (2*PHI_0)] = T * m² / Tm² = 1 --> without unit (phase)!
coeff = b_0 * a ** 2 / (2 * PHI_0) # Minus is gone because of negative z-direction
v_dim, u_dim = dim_uv
u = np.linspace(-(u_dim - 1), u_dim - 1, num=2 * u_dim - 1)
v = np.linspace(-(v_dim - 1), v_dim - 1, num=2 * v_dim - 1)
uu, vv = np.meshgrid(u, v)
self.u = fft.empty(self.dim_kern, dtype=fft.FLOAT)
self.v = fft.empty(self.dim_kern, dtype=fft.FLOAT)
self.u[...] = coeff * self._get_elementary_phase(geometry, uu, vv, a)
self.v[...] = coeff * -self._get_elementary_phase(geometry, vv, uu, a)
# Include perturbed reference wave:
if prw_vec is not None:
uu += prw_vec[1]
vv += prw_vec[0]
self.u[...] -= coeff * self._get_elementary_phase(geometry, uu, vv, a)
self.v[...] -= coeff * -self._get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
self.u_fft = fft.rfftn(self.u, self.dim_pad)
self.v_fft = fft.rfftn(self.v, self.dim_pad)
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, geometry=%r)' % \
(self.__class__, self.a, self.dim_uv, self.geometry)
def __str__(self):
self._log.debug('Calling __str__')
return 'Kernel(a=%s, dim_uv=%s, geometry=%s)' % \
(self.a, self.dim_uv, self.geometry)
def _get_elementary_phase(self, geometry, n, m, a):
self._log.debug('Calling _get_elementary_phase')
if geometry == 'disc':
in_or_out = ~ np.logical_and(n == 0, m == 0)
return m / (n ** 2 + m ** 2 + 1E-30) * in_or_out
elif geometry == 'slab':
def _F_a(n, m):
A = np.log(a ** 2 * (n ** 2 + m ** 2))
B = np.arctan(n / m)
return n * A - 2 * n + 2 * m * B
return 0.5 * (_F_a(m - 0.5, n - 0.5) - _F_a(m + 0.5, n - 0.5) -
_F_a(m - 0.5, n + 0.5) + _F_a(m + 0.5, n + 0.5))
def print_info(self):
"""Print information about the kernel.
Returns
-------
None
"""
self._log.debug('Calling log_info')
print('Shape of the FOV :', self.dim_uv)
print('Shape of the Kernel:', self.dim_kern)
print('Zero-padded shape :', self.dim_pad)
print('Shape of the FFT :', self.dim_fft)
print('Slice for the phase:', self.slice_phase)
print('Slice for the magn.:', self.slice_mag)
print('Grid spacing: {} nm'.format(self.a))
print('Geometry:', self.geometry)
pyramid/magcreator/__init__.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing functionality for creating magnetic distributions."""
from . import shapes
from . import examples
from .magcreator import *
__all__ = ['shapes', 'examples']
__all__.extend(magcreator.__all__)
pyramid/magcreator/examples.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Provide simple examples for magnetic distributions."""
import logging
import numpy as np
import random as rnd
from . import magcreator as mc
from . import shapes
from ..fielddata import VectorData
__all__ = ['pyramid_logo', 'singularity', 'homog_pixel', 'homog_slab', 'homog_disc',
'homog_sphere', 'homog_filament', 'homog_alternating_filament',
'homog_array_sphere_disc_slab', 'homog_random_pixels', 'homog_random_slabs',
'vortex_slab', 'vortex_disc', 'vortex_alternating_discs', 'vortex_sphere',
'vortex_horseshoe', 'core_shell_disc', 'core_shell_sphere']
_log = logging.getLogger(__name__)
def pyramid_logo(a=1., dim=(1, 256, 256), phi=np.pi / 2, theta=np.pi / 2):
"""Create pyramid logo."""
_log.debug('Calling pyramid_logo')
mag_shape = np.zeros(dim)
x = range(dim[2])
y = range(dim[1])
xx, yy = np.meshgrid(x, y)
bottom = (yy >= 0.25 * dim[1])
left = (yy <= 0.75 / 0.5 * dim[1] / dim[2] * xx)
right = np.fliplr(left)
mag_shape[0, ...] = np.logical_and(np.logical_and(left, right), bottom)
return VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
def singularity(a=1., dim=(5, 5, 5), center=None):
"""Create magnetic singularity."""
_log.debug('Calling singularity')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
zz, yy, xx = np.indices(dim)
magnitude = np.array((xx - center[2], yy - center[1], zz - center[0])).astype(float)
# magnitude /= np.sqrt((magnitude ** 2 + 1E-30).sum(axis=0))
return VectorData(a, magnitude)
def homog_pixel(a=1., dim=(1, 9, 9), pixel=None, phi=np.pi / 4, theta=np.pi / 2):
"""Create single magnetised slab."""
_log.debug('Calling homog_pixel')
if pixel is None:
pixel = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
mag_shape = shapes.pixel(dim, pixel)
return VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
def homog_slab(a=1., dim=(32, 32, 32), center=None, width=None, phi=np.pi / 4, theta=np.pi / 4):
"""Create homogeneous slab magnetisation distribution."""
_log.debug('Calling homog_slab')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if width is None:
width = (dim[0] // 8, dim[1] // 2, dim[2] // 4)
mag_shape = shapes.slab(dim, center, width)
return VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
def homog_disc(a=1., dim=(32, 32, 32), center=None, radius=None, height=None,
phi=np.pi / 4, theta=np.pi / 4):
"""Create homogeneous disc magnetisation distribution."""
_log.debug('Calling homog_disc')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius is None:
radius = dim[2] // 4
if height is None:
height = dim[0] // 2
mag_shape = shapes.disc(dim, center, radius, height)
return VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
def homog_sphere(a=1., dim=(32, 32, 32), center=None, radius=None, phi=np.pi / 4, theta=np.pi / 4):
"""Create homogeneous sphere magnetisation distribution."""
_log.debug('Calling homog_sphere')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius is None:
radius = dim[2] // 4
mag_shape = shapes.sphere(dim, center, radius)
return VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
def homog_filament(a=1., dim=(1, 21, 21), pos=None, phi=np.pi / 2, theta=np.pi / 2):
"""Create magnetisation distribution of a single magnetised filaments."""
_log.debug('Calling homog_filament')
if pos is None:
pos = (0, dim[1] // 2)
mag_shape = shapes.filament(dim, pos)
return VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
def homog_alternating_filament(a=1., dim=(1, 21, 21), spacing=5, phi=np.pi / 2, theta=np.pi / 2):
"""Create magnetisation distribution of alternating filaments."""
_log.debug('Calling homog_alternating_filament')
count = int((dim[1] - 1) / spacing) + 1
magdata = VectorData(a, np.zeros((3,) + dim))
for i in range(count):
pos = i * spacing
mag_shape = shapes.filament(dim, (0, pos))
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
phi *= -1 # Switch the angle
return magdata
def homog_array_sphere_disc_slab(a=1., dim=(64, 128, 128), center_sp=(32, 96, 64), radius_sp=24,
center_di=(32, 32, 96), radius_di=24, height_di=24,
center_sl=(32, 32, 32), width_sl=(48, 48, 48)):
"""Create array of several magnetisation distribution (sphere, disc and slab)."""
_log.debug('Calling homog_array_sphere_disc_slab')
mag_shape_sphere = shapes.sphere(dim, center_sp, radius_sp)
mag_shape_disc = shapes.disc(dim, center_di, radius_di, height_di)
mag_shape_slab = shapes.slab(dim, center_sl, width_sl)
magdata = VectorData(a, mc.create_mag_dist_homog(mag_shape_sphere, np.pi))
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape_disc, np.pi / 2))
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape_slab, np.pi / 4))
return magdata
def homog_random_pixels(a=1., dim=(1, 64, 64), count=10, rnd_seed=24):
"""Create random magnetised pixels."""
_log.debug('Calling homog_random_pixels')
rnd.seed(rnd_seed)
magdata = VectorData(a, np.zeros((3,) + dim))
for i in range(count):
pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
mag_shape = shapes.pixel(dim, pixel)
phi = 2 * np.pi * rnd.random()
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape, phi))
return magdata
def homog_random_slabs(a=1., dim=(1, 64, 64), count=10, width_max=5, rnd_seed=2):
"""Create random magnetised slabs."""
_log.debug('Create homog_random_slabs')
rnd.seed(rnd_seed)
magdata = VectorData(a, np.zeros((3,) + dim))
for i in range(count):
width = (1, rnd.randint(1, width_max), rnd.randint(1, width_max))
center = (rnd.randrange(int(width[0] / 2), dim[0] - int(width[0] / 2)),
rnd.randrange(int(width[1] / 2), dim[1] - int(width[1] / 2)),
rnd.randrange(int(width[2] / 2), dim[2] - int(width[2] / 2)))
mag_shape = shapes.slab(dim, center, width)
phi = 2 * np.pi * rnd.random()
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape, phi))
return magdata
def vortex_slab(a=1., dim=(32, 32, 32), center=None, width=None, axis='z'):
"""Create vortex slab magnetisation distribution."""
_log.debug('Calling vortex_slab')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if width is None:
width = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
mag_shape = shapes.slab(dim, center, width)
magnitude = mc.create_mag_dist_vortex(mag_shape, center, axis)
return VectorData(a, magnitude)
def vortex_disc(a=1., dim=(32, 32, 32), center=None, radius=None, height=None, axis='z'):
"""Create vortex disc magnetisation distribution."""
_log.debug('Calling vortex_disc')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius is None:
radius = dim[2] // 4
if height is None:
height = dim[0] // 2
mag_shape = shapes.disc(dim, center, radius, height, axis)
magnitude = mc.create_mag_dist_vortex(mag_shape, center, axis)
return VectorData(a, magnitude)
def vortex_alternating_discs(a=1., dim=(80, 32, 32), count=8):
"""Create pillar of alternating vortex disc magnetisation distributions."""
_log.debug('Calling vortex_alternating_discs')
segment_height = dim[0] // (count + 2)
magdata = VectorData(a, np.zeros((3,) + dim))
for i in range(count):
axis = 'z' if i % 2 == 0 else '-z'
center = (segment_height * (i + 1 + 0.5), dim[1] // 2, dim[2] // 2)
radius = dim[2] // 4
height = segment_height
mag_shape = shapes.disc(dim, center=center, radius=radius, height=height)
mag_amp = mc.create_mag_dist_vortex(mag_shape=mag_shape, center=center, axis=axis)
magdata += VectorData(1, mag_amp)
return magdata
def vortex_sphere(a=1., dim=(32, 32, 32), center=None, radius=None, axis = 'z'):
"""Create vortex sphere magnetisation distribution."""
_log.debug('Calling vortex_sphere')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius is None:
radius = dim[2] // 4
mag_shape = shapes.sphere(dim, center, radius)
magnitude = mc.create_mag_dist_vortex(mag_shape, center, axis)
return VectorData(a, magnitude)
def vortex_horseshoe(a=1., dim=(16, 64, 64), center=None, radius_core=None,
radius_shell=None, height=None):
"""Create magnetic horseshoe vortex magnetisation distribution."""
_log.debug('Calling vortex_horseshoe')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius_core is None:
radius_core = dim[1] // 8
if radius_shell is None:
radius_shell = dim[1] // 4
if height is None:
height = dim[0] // 2
mag_shape_core = shapes.disc(dim, center, radius_core, height)
mag_shape_outer = shapes.disc(dim, center, radius_shell, height)
mag_shape_horseshoe = np.logical_xor(mag_shape_outer, mag_shape_core)
mag_shape_horseshoe[:, dim[1] // 2:, :] = False
return VectorData(a, mc.create_mag_dist_vortex(mag_shape_horseshoe))
def core_shell_disc(a=1., dim=(32, 32, 32), center=None, radius_core=None,
radius_shell=None, height=None, rate_core_to_shell=0.75):
"""Create magnetic core shell disc magnetisation distribution."""
_log.debug('Calling core_shell_disc')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius_core is None:
radius_core = dim[1] // 8
if radius_shell is None:
radius_shell = dim[1] // 4
if height is None:
height = dim[0] // 2
mag_shape_core = shapes.disc(dim, center, radius_core, height)
mag_shape_outer = shapes.disc(dim, center, radius_shell, height)
mag_shape_shell = np.logical_xor(mag_shape_outer, mag_shape_core)
magdata = VectorData(a, mc.create_mag_dist_vortex(mag_shape_shell)) * rate_core_to_shell
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
return magdata
def core_shell_sphere(a=1., dim=(32, 32, 32), center=None, radius_core=None,
radius_shell=None, rate_core_to_shell=0.75):
"""Create magnetic core shell sphere magnetisation distribution."""
_log.debug('Calling core_shell_sphere')
if center is None:
center = (dim[0] // 2, dim[1] // 2, dim[2] // 2)
if radius_core is None:
radius_core = dim[1] // 8
if radius_shell is None:
radius_shell = dim[1] // 4
mag_shape_sphere = shapes.sphere(dim, center, radius_shell)
mag_shape_disc = shapes.disc(dim, center, radius_core, height=dim[0])
mag_shape_core = np.logical_and(mag_shape_sphere, mag_shape_disc)
mag_shape_shell = np.logical_and(mag_shape_sphere, np.logical_not(mag_shape_core))
magdata = VectorData(a, mc.create_mag_dist_vortex(mag_shape_shell)) * rate_core_to_shell
magdata += VectorData(a, mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
return magdata
pyramid/magcreator/magcreator.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Create simple magnetic distributions.
The :mod:`~.magcreator` module is responsible for the creation of simple distributions of
magnetic moments. In the :mod:`~.shapes` module, you can find several general shapes for the
3-dimensional volume that should be magnetized (e.g. slabs, spheres, discs or single pixels).
These shapes are then used as input for the creating functions (or you could specify the
volume yourself as a 3-dimensional boolean matrix or a matrix with values in the range from 0 to 1,
which modifies the magnetization amplitude). The specified volume can either be magnetized
homogeneously with the :func:`~.create_mag_dist_homog` function by specifying the magnetization
direction, or in a vortex state with the :func:`~.create_mag_dist_vortex` by specifying the vortex
center.
"""
import logging
import numpy as np
from numpy import pi
__all__ = ['create_mag_dist_homog', 'create_mag_dist_vortex']
_log = logging.getLogger(__name__)
def create_mag_dist_homog(mag_shape, phi, theta=pi / 2):
"""Create a 3-dimensional magnetic distribution of a homogeneously magnetized object.
Parameters
----------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see :mod:`.~shapes` for examples).
phi : float
The azimuthal angle, describing the direction of the magnetized object.
theta : float, optional
The polar angle, describing the direction of the magnetized object.
The default is pi/2, which means, the z-component is zero.
Returns
-------
amplitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
"""
_log.debug('Calling create_mag_dist_homog')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
z_mag = np.ones(dim) * np.cos(theta) * mag_shape
y_mag = np.ones(dim) * np.sin(theta) * np.sin(phi) * mag_shape
x_mag = np.ones(dim) * np.sin(theta) * np.cos(phi) * mag_shape
return np.array([x_mag, y_mag, z_mag])
def create_mag_dist_vortex(mag_shape, center=None, axis='z'):
"""Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
----------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see :mod:`.~shapes`` for examples).
center : tuple (N=2 or N=3), optional
The vortex center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
is is discarded. Is set to the center of the field of view if not specified. The vortex
center has to be between two pixels.
axis : {'z', '-z', 'y', '-y', 'x', '-x'}, optional
The orientation of the vortex axis. The default is 'z'. Negative values invert the vortex
orientation.
Returns
-------
amplitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
Notes
-----
To avoid singularities, the vortex center should lie between the pixel centers (which
reside at coordinates with _.5 at the end), i.e. integer values should be used as center
coordinates (e.g. coordinate 1 lies between the first and the second pixel).
"""
_log.debug('Calling create_mag_dist_vortex')
dim = mag_shape.shape
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
assert center is None or len(center) in {2, 3}, \
'Vortex center has to be defined in 3D or 2D or not at all!'
if center is None:
center = (dim[1] / 2, dim[2] / 2)
sign = -1 if '-' in axis else 1
if axis in ('z', '-z'):
if len(center) == 3: # if a 3D-center is given, just take the x and y components
center = (center[1], center[2])
u = np.linspace(-center[1], dim[2] - 1 - center[1], dim[2]) + 0.5 # pixel center!
v = np.linspace(-center[0], dim[1] - 1 - center[0], dim[1]) + 0.5 # pixel center!
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu) - pi / 2, axis=0)
phi = np.tile(phi, (dim[0], 1, 1))
z_mag = np.zeros(dim)
y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * sign
x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * sign
elif axis in ('y', '-y'):
if len(center) == 3: # if a 3D-center is given, just take the x and z components
center = (center[0], center[2])
u = np.linspace(-center[1], dim[2] - 1 - center[1], dim[2]) + 0.5 # pixel center!
v = np.linspace(-center[0], dim[0] - 1 - center[0], dim[0]) + 0.5 # pixel center!
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu) - pi / 2, axis=1)
phi = np.tile(phi, (1, dim[1], 1))
z_mag = np.ones(dim) * np.sin(phi) * mag_shape * sign
y_mag = np.zeros(dim)
x_mag = np.ones(dim) * np.cos(phi) * mag_shape * sign
elif axis in ('x', '-x'):
if len(center) == 3: # if a 3D-center is given, just take the z and y components
center = (center[0], center[1])
u = np.linspace(-center[1], dim[1] - 1 - center[1], dim[1]) + 0.5 # pixel center!
v = np.linspace(-center[0], dim[0] - 1 - center[0], dim[0]) + 0.5 # pixel center!
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu) - pi / 2, axis=2)
phi = np.tile(phi, (1, 1, dim[2]))
z_mag = -np.ones(dim) * np.sin(phi) * mag_shape * sign
y_mag = -np.ones(dim) * np.cos(phi) * mag_shape * sign
x_mag = np.zeros(dim)
else:
raise ValueError('{} is not a valid argument (use x, -x, y, -y, z or -z)'.format(axis))
return np.array([x_mag, y_mag, z_mag])
pyramid/magcreator/shapes.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Provide simple shapes.
This module is a collection of some methods that return a 3-dimensional
matrix that represents the field volume and consists of boolean values.
This matrix is used in the functions of the :mod:`~.magcreator` module to create
:class:`~pyramid.fielddata.VectorData` objects which store the field information.
"""
import logging
import numpy as np
__all__ = ['slab', 'disc', 'ellipse', 'ellipsoid', 'sphere', 'filament', 'pixel']
_log = logging.getLogger(__name__)
def slab(dim, center, width):
"""Create the shape of a slab.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the slab in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the slab in pixel coordinates `(z, y, x)`.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling slab')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a tuple of length 3!'
zz, yy, xx = np.indices(dim) + 0.5
xx_shape = np.where(abs(xx - center[2]) <= width[2] / 2, True, False)
yy_shape = np.where(abs(yy - center[1]) <= width[1] / 2, True, False)
zz_shape = np.where(abs(zz - center[0]) <= width[0] / 2, True, False)
return np.logical_and(np.logical_and(xx_shape, yy_shape), zz_shape)
def disc(dim, center, radius, height, axis='z'):
"""Create the shape of a cylindrical disc in x-, y-, or z-direction.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
axis : {'z', 'y', 'x'}, optional
The orientation of the disc axis. The default is 'z'.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling disc')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
zz, yy, xx = np.indices(dim) + 0.5
xx -= center[2]
yy -= center[1]
zz -= center[0]
if axis == 'z':
uu, vv, ww = xx, yy, zz
elif axis == 'y':
uu, vv, ww = zz, xx, yy
elif axis == 'x':
uu, vv, ww = yy, zz, xx
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(axis))
return np.logical_and(np.where(np.hypot(uu, vv) <= radius, True, False),
np.where(abs(ww) <= height / 2, True, False))
def ellipse(dim, center, width, height, axis='z'):
"""Create the shape of an elliptical cylinder in x-, y-, or z-direction.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the ellipse in pixel coordinates `(z, y, x)`.
width : tuple (N=2)
Length of the two axes of the ellipse.
height : float
The height of the ellipse in pixel coordinates.
axis : {'z', 'y', 'x'}, optional
The orientation of the ellipse axis. The default is 'z'.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling ellipse')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (2,), 'Parameter width has to be a a tuple of length 2!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
zz, yy, xx = np.indices(dim) + 0.5
xx -= center[2]
yy -= center[1]
zz -= center[0]
if axis == 'z':
uu, vv, ww = xx, yy, zz
elif axis == 'y':
uu, vv, ww = xx, zz, yy
elif axis == 'x':
uu, vv, ww = yy, zz, xx
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(axis))
distance = np.hypot(uu / (width[1] / 2), vv / (width[0] / 2))
return np.logical_and(np.where(distance <= 1, True, False),
np.where(abs(ww) <= height / 2, True, False))
def ellipsoid(dim, center, width):
"""Create the shape of an ellipsoid.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the ellipsoid in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the ellipsoid `(z, y, x)`.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling ellipsoid')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a a tuple of length 3!'
zz, yy, xx = np.indices(dim) + 0.5
distance = np.sqrt(((xx - center[2]) / (width[2] / 2)) ** 2
+ ((yy - center[1]) / (width[1] / 2)) ** 2
+ ((zz - center[0]) / (width[0] / 2)) ** 2)
return np.where(distance <= 1, True, False)
def sphere(dim, center, radius):
"""Create the shape of a sphere.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the sphere in pixel coordinates `(z, y, x)`.
radius : float
The radius of the sphere in pixel coordinates.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling sphere')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
zz, yy, xx = np.indices(dim) + 0.5
distance = np.sqrt((xx - center[2]) ** 2 + (yy - center[1]) ** 2 + (zz - center[0]) ** 2)
return np.where(distance <= radius, True, False)
def filament(dim, pos, axis='y'):
"""Create the shape of a filament.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
pos : tuple (N=2)
The position of the filament in pixel coordinates `(coord1, coord2)`.
`coord1` and `coord2` stand for the two axes, which are perpendicular to `axis`. For
the default case (`axis = y`), it is `(coord1, coord2) = (z, x)`.
axis : {'y', 'x', 'z'}, optional
The orientation of the filament axis. The default is 'y'.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling filament')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pos) == (2,), 'Parameter pos has to be a tuple of length 2!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
shape = np.zeros(dim, dtype=bool)
if axis == 'z':
shape[:, pos[0], pos[1]] = True
elif axis == 'y':
shape[pos[0], :, pos[1]] = True
elif axis == 'x':
shape[pos[0], pos[1], :] = True
return shape
def pixel(dim, pixel):
"""Create the shape of a single pixel.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
pixel : tuple (N=3)
The coordinates of the pixel `(z, y, x)`.
Returns
-------
shape : :class:`~numpy.ndarray` (N=3)
The shape as a 3D-array.
"""
_log.debug('Calling pixel')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pixel) == (3,), 'Parameter pixel has to be a tuple of length 3!'
shape = np.zeros(dim, dtype=bool)
shape[pixel] = True
return shape
pyramid/phasemap.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.PhaseMap` class for storing phase map data."""
import logging
from numbers import Number
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.ticker import MaxNLocator, FuncFormatter
from mpl_toolkits.mplot3d import Axes3D
from scipy.ndimage.interpolation import zoom
from . import colors
__all__ = ['PhaseMap']
class PhaseMap(object):
"""Class for storing phase map data.
Represents 2-dimensional phase maps. The phase information itself is stored as a 2-dimensional
matrix in `phase`, but can also be accessed as a vector via `phase_vec`. :class:`~.PhaseMap`
objects support negation, arithmetic operators (``+``, ``-``, ``*``) and their augmented
counterparts (``+=``, ``-=``, ``*=``), with numbers and other :class:`~.PhaseMap`
objects, if their dimensions and grid spacings match. It is possible to load data from HDF5
or textfiles or to save the data in these formats. Methods for plotting the phase or a
corresponding holographic contour map are provided. Holographic contour maps are created by
taking the cosine of the (optionally amplified) phase and encoding the direction of the
2-dimensional gradient via color. The directional encoding can be seen by using the
:func:`~.make_color_wheel` function. Use the :func:`~.plot_combined` function to plot the
phase map and the holographic contour map next to each other.
Attributes
----------
a: float
The grid spacing in nm.
phase: :class:`~numpy.ndarray` (N=2)
Array containing the phase shift.
mask: :class:`~numpy.ndarray` (boolean, N=2, optional)
Mask which determines the projected magnetization distribution, gotten from MIP images or
otherwise acquired. Defaults to an array of ones (all pixels are considered).
confidence: :class:`~numpy.ndarray` (N=2, optional)
Confidence array which determines the trust of specific regions of the phasemap. A value
of 1 means the pixel is trustworthy, a value of 0 means it is not. Defaults to an array of
ones (full trust for all pixels). Can be used for the construction of Se_inv.
"""
_log = logging.getLogger(__name__)
UNITDICT = {u'rad': 1E0,
u'mrad': 1E3,
u'µrad': 1E6}
CDICT = {'red': [(0.00, 1.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 1.0, 1.0),
(0.75, 0.0, 0.0),
(1.00, 0.0, 1.0)],
'green': [(0.00, 0.0, 0.0),
(0.25, 0.0, 0.0),
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 1.0)],
'blue': [(0.00, 1.0, 1.0),
(0.25, 0.0, 0.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 1.0, 1.0)]}
CDICT_INV = {'red': [(0.00, 0.0, 1.0),
(0.25, 0.0, 0.0),
(0.50, 0.0, 0.0),
(0.75, 1.0, 1.0),
(1.00, 1.0, 0.0)],
'green': [(0.00, 1.0, 1.0),
(0.25, 1.0, 1.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 1.0, 0.0)],
'blue': [(0.00, 0.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 0.0)]}
HOLO_CMAP = LinearSegmentedColormap('my_colormap', CDICT, 256)
HOLO_CMAP_INV = LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
@property
def a(self):
"""Grid spacing in nm."""
return self._a
@a.setter
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = float(a)
@property
def dim_uv(self):
"""Dimensions of the grid."""
return self._dim_uv
@property
def phase(self):
"""Array containing the phase shift."""
return self._phase
@phase.setter
def phase(self, phase):
assert isinstance(phase, np.ndarray), 'Phase has to be a numpy array!'
assert len(phase.shape) == 2, 'Phase has to be 2-dimensional!'
self._phase = phase.astype(dtype=np.float32)
self._dim_uv = phase.shape
@property
def phase_vec(self):
"""Vector containing the phase shift."""
return np.reshape(self.phase, -1)
@phase_vec.setter
def phase_vec(self, phase_vec):
assert isinstance(phase_vec, np.ndarray), 'Vector has to be a numpy array!'
assert np.size(phase_vec) == np.prod(self.dim_uv), 'Vector size has to match phase!'
self.phase = phase_vec.reshape(self.dim_uv)
@property
def mask(self):
"""Mask which determines the projected magnetization distribution"""
return self._mask
@mask.setter
def mask(self, mask):
if mask is not None:
assert mask.shape == self.phase.shape, 'Mask and phase dimensions must match!!'
else:
mask = np.ones_like(self.phase, dtype=bool)
self._mask = mask.astype(np.bool)
@property
def confidence(self):
"""Confidence array which determines the trust of specific regions of the phasemap."""
return self._confidence
@confidence.setter
def confidence(self, confidence):
if confidence is not None:
assert confidence.shape == self.phase.shape, \
'Confidence and phase dimensions must match!'
else:
confidence = np.ones_like(self.phase)
self._confidence = confidence.astype(dtype=np.float32)
def __init__(self, a, phase, mask=None, confidence=None):
self._log.debug('Calling __init__')
self.a = a
self.phase = phase
self.mask = mask
self.confidence = confidence
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, phase=%r, mask=%r, confidence=%r)' % \
(self.__class__, self.a, self.phase, self.mask, self.confidence)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMap(a=%s, dim_uv=%s, mask=%s)' % (self.a, self.dim_uv, not np.all(self.mask))
def __neg__(self): # -self
self._log.debug('Calling __neg__')
return PhaseMap(self.a, -self.phase, self.mask, self.confidence)
def __add__(self, other): # self + other
self._log.debug('Calling __add__')
assert isinstance(other, (PhaseMap, Number)), \
'Only PhaseMap objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, PhaseMap):
self._log.debug('Adding two PhaseMap objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.phase.shape == self.dim_uv, \
'Added field has to have the same dimensions!'
mask_comb = np.logical_or(self.mask, other.mask) # masks combine
conf_comb = (self.confidence + other.confidence) / 2 # confidence averaged!
return PhaseMap(self.a, self.phase + other.phase, mask_comb, conf_comb)
else: # other is a Number
self._log.debug('Adding an offset')
return PhaseMap(self.a, self.phase + other, self.mask, self.confidence)
def __sub__(self, other): # self - other
self._log.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
self._log.debug('Calling __mul__')
assert (isinstance(other, Number) or
(isinstance(other, np.ndarray) and other.shape == self.dim_uv)), \
'PhaseMap objects can only be multiplied by scalar numbers or fitting arrays!'
return PhaseMap(self.a, self.phase * other, self.mask, self.confidence)
def __truediv__(self, other): # self / other
self._log.debug('Calling __truediv__')
assert (isinstance(other, Number) or
(isinstance(other, np.ndarray) and other.shape == self.dim_uv)), \
'PhaseMap objects can only be divided by scalar numbers or fitting arrays!'
return PhaseMap(self.a, self.phase / other, self.mask, self.confidence)
def __floordiv__(self, other): # self // other
self._log.debug('Calling __floordiv__')
assert (isinstance(other, Number) or
(isinstance(other, np.ndarray) and other.shape == self.dim_uv)), \
'PhaseMap objects can only be divided by scalar numbers or fitting arrays!'
return PhaseMap(self.a, self.phase // other, self.mask, self.confidence)
def __radd__(self, other): # other + self
self._log.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
self._log.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
self._log.debug('Calling __rmul__')
return self.__mul__(other)
def __iadd__(self, other): # self += other
self._log.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
self._log.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
self._log.debug('Calling __imul__')
return self.__mul__(other)
def __itruediv__(self, other): # self /= other
self._log.debug('Calling __itruediv__')
return self.__truediv__(other)
def __ifloordiv__(self, other): # self //= other
self._log.debug('Calling __ifloordiv__')
return self.__floordiv__(other)
def __array__(self, dtype=None):
if dtype:
return self.phase.astype(dtype)
else:
return self.phase
def __array_wrap__(self, array, _=None): # _ catches the context, which is not used.
return PhaseMap(self.a, array, self.mask, self.confidence)
def copy(self):
"""Returns a copy of the :class:`~.PhaseMap` object
Returns
-------
phasemap: :class:`~.PhaseMap`
A copy of the :class:`~.PhaseMap`.
"""
self._log.debug('Calling copy')
return PhaseMap(self.a, self.phase.copy(), self.mask.copy(),
self.confidence.copy())
def scale_down(self, n=1):
"""Scale down the phase map by averaging over two pixels along each axis.
Parameters
----------
n : int, optional
Number of times the phase map is scaled down. The default is 1.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
Only possible, if each axis length is a power of 2!
"""
self._log.debug('Calling scale_down')
assert n > 0 and isinstance(n, int), 'n must be a positive integer!'
self.a *= 2 ** n
for t in range(n):
# Pad if necessary:
pv, pu = (0, self.dim_uv[0] % 2), (0, self.dim_uv[1] % 2)
if pv != 0 or pu != 0:
self.pad((pv, pu), mode='edge')
# Create coarser grid for the magnetization:
dim_uv = self.dim_uv
self.phase = self.phase.reshape((dim_uv[0] // 2, 2,
dim_uv[1] // 2, 2)).mean(axis=(3, 1))
mask = self.mask.reshape(dim_uv[0] // 2, 2, dim_uv[1] // 2, 2)
self.mask = mask[:, 0, :, 0] & mask[:, 1, :, 0] & mask[:, 0, :, 1] & mask[:, 1, :, 1]
self.confidence = self.confidence.reshape(dim_uv[0] // 2, 2,
dim_uv[1] // 2, 2).mean(axis=(3, 1))
def scale_up(self, n=1, order=0):
"""Scale up the phase map using spline interpolation of the requested order.
Parameters
----------
n : int, optional
Power of 2 with which the grid is scaled. Default is 1, which means every axis is
increased by a factor of ``2**1 = 2``.
order : int, optional
The order of the spline interpolation, which has to be in the range between 0 and 5
and defaults to 0.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
"""
self._log.debug('Calling scale_up')
assert n > 0 and isinstance(n, int), 'n must be a positive integer!'
assert 5 > order >= 0 and isinstance(order, int), \
'order must be a positive integer between 0 and 5!'
self.a /= 2 ** n
self.phase = zoom(self.phase, zoom=2 ** n, order=order)
self.mask = zoom(self.mask, zoom=2 ** n, order=0)
self.confidence = zoom(self.confidence, zoom=2 ** n, order=order)
def pad(self, pad_values, mode='constant', masked=False, **kwds):
"""Pad the current phase map with zeros for each individual axis.
Parameters
----------
pad_values : tuple of int
Number of zeros which should be padded. Provided as a tuple where each entry
corresponds to an axis. An entry can be one int (same padding for both sides) or again
a tuple which specifies the pad values for both sides of the corresponding axis.
mode: string or function
A string values or a user supplied function. ‘constant’ pads with zeros. ‘edge’ pads
with the edge values of array. See the numpy pad function for an in depth guide.
masked: boolean, optional
Determines if the padded areas should be masked or not. `True` creates a 'buffer
zone' for the magnetization distribution in the reconstruction. Default is `False`
Returns
-------
None
Notes
-----
Acts in place and changes dimensions accordingly.
The confidence of the padded areas is set to zero!
"""
self._log.debug('Calling pad')
assert len(pad_values) == 2, 'Pad values for each dimension have to be provided!'
pval = np.zeros(4, dtype=np.int)
for i, values in enumerate(pad_values):
assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
pval[2 * i:2 * (i + 1)] = values
self.phase = np.pad(self.phase, ((pval[0], pval[1]),
(pval[2], pval[3])), mode=mode, **kwds)
self.confidence = np.pad(self.confidence, ((pval[0], pval[1]),
(pval[2], pval[3])), mode=mode, **kwds)
if masked:
mask_kwds = {'mode': 'constant', 'constant_values': True}
else:
mask_kwds = {'mode': mode}
self.mask = np.pad(self.mask, ((pval[0], pval[1]), (pval[2], pval[3])), **mask_kwds)
def crop(self, crop_values):
"""Pad the current phase map with zeros for each individual axis.
Parameters
----------
crop_values : tuple of int
Number of zeros which should be cropped. Provided as a tuple where each entry
corresponds to an axis. An entry can be one int (same cropping for both sides) or again
a tuple which specifies the crop values for both sides of the corresponding axis.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions accordingly.
"""
self._log.debug('Calling crop')
assert len(crop_values) == 2, 'Crop values for each dimension have to be provided!'
cv = np.zeros(4, dtype=np.int)
for i, values in enumerate(crop_values):
assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
cv[2 * i:2 * (i + 1)] = values
cv *= np.resize([1, -1], len(cv))
cv = np.where(cv == 0, None, cv)
self.phase = self.phase[cv[0]:cv[1], cv[2]:cv[3]]
self.mask = self.mask[cv[0]:cv[1], cv[2]:cv[3]]
self.confidence = self.confidence[cv[0]:cv[1], cv[2]:cv[3]]
@classmethod
def from_signal(cls, signal):
"""Convert a :class:`~hyperspy.signals.Image` object to a :class:`~.PhaseMap` object.
Parameters
----------
signal: :class:`~hyperspy.signals.Image`
The :class:`~hyperspy.signals.Image` object which should be converted to a PhaseMap.
Returns
-------
phasemap: :class:`~.PhaseMap`
A :class:`~.PhaseMap` object containing the loaded data.
Notes
-----
This method recquires the hyperspy package!
"""
cls._log.debug('Calling from_signal')
# Extract phase:
phase = signal.data
# Extract properties:
a = signal.axes_manager.signal_axes[0].scale
try:
mask = signal.metadata.Signal.mask
confidence = signal.metadata.Signal.confidence
except AttributeError:
mask = None
confidence = None
return cls(a, phase, mask, confidence)
def to_signal(self):
"""Convert :class:`~.PhaseMap` data into a HyperSpy Signal.
Returns
-------
signal: :class:`~hyperspy.signals.Signal2D`
Representation of the :class:`~.PhaseMap` object as a HyperSpy Signal.
Notes
-----
This method recquires the hyperspy package!
"""
self._log.debug('Calling to_signal')
try: # Try importing HyperSpy:
import hyperspy.api as hs
except ImportError:
self._log.error('This method recquires the hyperspy package!')
return
# Create signal:
signal = hs.signals.Signal2D(self.phase)
# Set axes:
signal.axes_manager.signal_axes[0].name = 'x-axis'
signal.axes_manager.signal_axes[0].units = 'nm'
signal.axes_manager.signal_axes[0].scale = self.a
signal.axes_manager.signal_axes[1].name = 'y-axis'
signal.axes_manager.signal_axes[1].units = 'nm'
signal.axes_manager.signal_axes[1].scale = self.a
# Set metadata:
signal.metadata.Signal.title = 'PhaseMap'
signal.metadata.Signal.unit = 'rad'
signal.metadata.Signal.mask = self.mask
signal.metadata.Signal.confidence = self.confidence
# Create and return EMD:
return signal
def save(self, filename, save_mask=False, save_conf=False, pyramid_format=True, **kwargs):
"""Saves the phasemap in the specified format.
The function gets the format from the extension:
- hdf5 for HDF5.
- rpl for Ripple (useful to export to Digital Micrograph).
- unf for SEMPER unf binary format.
- txt format.
- Many image formats such as png, tiff, jpeg...
If no extension is provided, 'hdf5' is used. Most formats are
saved with the HyperSpy package (internally the phasemap is first
converted to a HyperSpy Signal.
Each format accepts a different set of parameters. For details
see the specific format documentation.
Parameters
----------
filename: str, optional
Name of the file which the phasemap is saved into. The extension
determines the saving procedure.
save_mask: boolean, optional
If True, the `mask` is saved, too. For all formats, except HDF5, a separate file will
be created. HDF5 always saves the `mask` in the metadata, independent of this flag. The
default is False.
save_conf: boolean, optional
If True, the `confidence` is saved, too. For all formats, except HDF5, a separate file
will be created. HDF5 always saves the `confidence` in the metadata, independent of
this flag. The default is False
pyramid_format: boolean, optional
Only used for saving to '.txt' files. If this is True, the grid spacing is saved
in an appropriate header. Otherwise just the phase is written with the
corresponding `kwargs`.
"""
from .file_io.io_phasemap import save_phasemap
save_phasemap(self, filename, save_mask, save_conf, pyramid_format, **kwargs)
def plot_phase(self, title='Phase Map', cbar_title=None, unit='rad', cmap='RdBu', limit=None,
norm=None, axis=None, cbar=True, show_mask=True, show_conf=True,
sigma_clip=None, interpolation='none'):
"""Display the phasemap as a colormesh.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Phase Map'.
cbar_title : string, optional
The title of the colorbar.
unit: {'rad', 'mrad', 'µrad'}, optional
The plotting unit of the phase map. The phase is scaled accordingly before plotting.
always in `rad`.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
cbar : bool, optional
A switch determining if the colorbar should be plotted or not. Default is True.
show_mask : bool, optional
A switch determining if the mask should be plotted or not. Default is True.
show_conf : float, optional
A switch determining if the confidence should be plotted or not. Default is True.
sigma_clip : int, optional
If this is not `None`, the values outside `sigma_clip` times the standard deviation
will be clipped for the calculation of the plotting `limit`.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
Returns
-------
axis, cbar: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted and the colorbar.
"""
self._log.debug('Calling plot_phase')
# Take units into consideration:
phase = self.phase * self.UNITDICT[unit]
# Calculate limit if necessary:
if limit is None:
phase_l = phase
# Clip non-trustworthy regions for the limit calculation:
if show_conf:
phase_trust = np.where(self.confidence > 0.9, phase_l, np.nan)
phase_min, phase_max = np.nanmin(phase_trust), np.nanmax(phase_trust)
phase_l = np.clip(phase_l, phase_min, phase_max)
# Cut outlier beyond a certain sigma-margin:
if sigma_clip is not None:
outlier = np.abs(phase_l - np.mean(phase_l)) < sigma_clip * np.std(phase_l)
phase_sigma = np.where(outlier, phase_l, np.nan)
phase_min, phase_max = np.nanmin(phase_sigma), np.nanmax(phase_sigma)
phase_l = np.clip(phase_l, phase_min, phase_max)
# Calculate the limit:
limit = np.max(np.abs(phase_l))
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure(figsize=(7, 7))
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the phasemap:
im = axis.imshow(phase, cmap=cmap, vmin=-limit, vmax=limit, interpolation=interpolation,
norm=norm, origin='lower', extent=(0, self.dim_uv[1], 0, self.dim_uv[0]))
if show_mask or show_conf:
vv, uu = np.indices(self.dim_uv) + 0.5
if show_conf and not np.all(self.confidence == 1.0):
colormap = colors.transparent_cmap
axis.imshow(self.confidence, cmap=colormap, interpolation=interpolation,
origin='lower', extent=(0, self.dim_uv[1], 0, self.dim_uv[0]))
if show_mask and not np.all(self.mask): # Plot mask if desired and not trivial!
axis.contour(uu, vv, self.mask, levels=[0.5], colors='k', linestyles='dotted',
linewidths=2)
# Set the axes ticks and labels:
if self.dim_uv[0] >= self.dim_uv[1]:
u_bin, v_bin = np.max((2, np.floor(9 * self.dim_uv[1] / self.dim_uv[0]))), 9
else:
u_bin, v_bin = 9, np.max((2, np.floor(9 * self.dim_uv[0] / self.dim_uv[1])))
axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x * self.a)))
axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x * self.a)))
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlim(0, self.dim_uv[1])
axis.set_ylim(0, self.dim_uv[0])
axis.set_xlabel('u-axis [nm]', fontsize=15)
axis.set_ylabel('v-axis [nm]', fontsize=15)
# # Add colorbar:
if cbar:
fig = plt.gcf()
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.82, 0.15, 0.02, 0.7])
cbar = fig.colorbar(im, cax=cbar_ax)
cbar.ax.tick_params(labelsize=14)
if cbar_title is None:
cbar_title = u'phase shift [{}]'.format(unit)
cbar.set_label(cbar_title, fontsize=15)
# Return plotting axis:
return axis
def plot_phase3d(self, title='Phase Map', unit='rad', cmap='RdBu'):
"""Display the phasemap as a 3D surface with contourplots.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Phase Map'.
unit: {'rad', 'mrad', 'µrad'}, optional
The plotting unit of the phase map. The phase is scaled accordingly before plotting.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
"""
self._log.debug('Calling plot_phase3d')
# Take units into consideration:
phase = self.phase * self.UNITDICT[unit]
# Create figure and axis:
fig = plt.figure()
axis = Axes3D(fig)
# Plot surface and contours:
vv, uu = np.indices(self.dim_uv)
axis.plot_surface(uu, vv, phase, rstride=4, cstride=4, alpha=0.7, cmap=cmap,
linewidth=0, antialiased=False)
axis.contourf(uu, vv, phase, 15, zdir='z', offset=np.min(phase), cmap=cmap)
axis.set_title(title)
axis.view_init(45, -135)
axis.set_xlabel('u-axis [px]')
axis.set_ylabel('v-axis [px]')
axis.set_zlabel('phase shift [{}]'.format(unit))
# Return plotting axis:
return axis
def plot_holo(self, title=None, gain='auto', axis=None, hue_mode='triadic',
interpolation='none'):
"""Display the color coded holography image.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Contour Map (gain: %g)' % gain.
gain : float or 'auto', optional
The gain factor for determining the number of contour lines. The default is 'auto',
which means that the gain will be determined automatically to look pretty.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method. No interpolation is used in the default case.
hue_mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme. Use either a triadic or tetradic
scheme (see the according colormaps for more information).
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
"""
self._log.debug('Calling plot_holo')
# Calculate gain if 'auto' is selected:
if gain == 'auto':
gain = 4 * 2 * np.pi / (np.abs(self.phase).max() + 1E-30)
# Set title if not set:
if title is None:
title = 'Holographic contour Map (gain: {:g})'.format(gain)
# Calculate the holography image intensity:
holo = np.cos(gain * self.phase)
holo += 1 # Shift to positive values
holo /= 2 # Rescale to [0, 1]
# Calculate the phase gradients and calculate colors:
# B = rot(A) --> B_x = grad_y(A_z), B_y = -grad_x(A_z); phi_m ~ -int(A_z)
# sign switch --> B_x = -grad_y(phi_m), B_y = grad_x(phi_m)
grad_x, grad_y = np.gradient(self.phase, self.a, self.a)
rgb = colors.rgb_from_vector(grad_x, -grad_y, np.zeros_like(grad_x), mode=hue_mode)
rgb = (holo.T * rgb.T).T.astype(np.uint8)
holo_image = Image.fromarray(rgb)
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the image and set axes:
axis.imshow(holo_image, origin='lower', interpolation=interpolation,
extent=(0, self.dim_uv[1], 0, self.dim_uv[0]))
# Set the title and the axes labels:
axis.set_title(title)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('u-axis [px]', fontsize=15)
axis.set_ylabel('v-axis [px]', fontsize=15)
if self.dim_uv[0] >= self.dim_uv[1]:
u_bin, v_bin = np.max((2, np.floor(9 * self.dim_uv[1] / self.dim_uv[0]))), 9
else:
u_bin, v_bin = 9, np.max((2, np.floor(9 * self.dim_uv[0] / self.dim_uv[1])))
axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
# Return plotting axis:
return axis
def plot_combined(self, sup_title='Combined Plot', phase_title='Phase Map', holo_title=None,
cbar_title=None, unit='rad', cmap='RdBu', limit=None, norm=None, gain='auto',
interpolation='none', cbar=True, show_mask=True, show_conf=True,
hue_mode='triadic'):
"""Display the phase map and the resulting color coded holography image in one plot.
Parameters
----------
sup_title : string, optional
The super title of the plot. The default is 'Combined Plot'.
phase_title : string, optional
The title of the phase map.
holo_title : string, optional
The title of the holographic contour map
cbar_title : string, optional
The title of the colorbar.
unit: {'rad', 'mrad', 'µrad'}, optional
The plotting unit of the phase map. The phase is scaled accordingly before plotting.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
gain : float or 'auto', optional
The gain factor for determining the number of contour lines. The default is 'auto',
which means that the gain will be determined automatically to look pretty.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
cbar : bool, optional
A switch determining if the colorbar should be plotted or not. Default is True.
show_mask : bool, optional
A switch determining if the mask should be plotted or not. Default is True.
show_conf : float, optional
A switch determining if the confidence should be plotted or not. Default is True.
hue_mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme. Use either a triadic or tetradic
scheme (see the according colormaps for more information).
Returns
-------
phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
The axes on which the graphs are plotted.
"""
self._log.debug('Calling plot_combined')
# Create combined plot and set title:
fig = plt.figure(figsize=(15, 7))
fig.suptitle(sup_title, fontsize=20)
# Plot holography image:
holo_axis = fig.add_subplot(1, 2, 1, aspect='equal')
self.plot_holo(title=holo_title, gain=gain, axis=holo_axis, interpolation=interpolation,
hue_mode=hue_mode)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
self.plot_phase(title=phase_title, cbar_title=cbar_title, unit=unit, cmap=cmap,
limit=limit, norm=norm, axis=phase_axis, cbar=cbar,
show_mask=show_mask, show_conf=show_conf)
# Return the plotting axes:
return phase_axis, holo_axis
pyramid/phasemapper.py deleted 100644 → 0
View file @ 68453aaa
# -*- 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
from . import fft
from .fielddata import VectorData, ScalarData
from .phasemap import PhaseMap
__all__ = ['PhaseMapperRDFC', 'PhaseMapperFDFC', 'PhaseMapperMIP', 'PhaseMapperCharge']
_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
EPS_0 = 8.8542E-12 # electrical field constant
class PhaseMapper(object, metaclass=abc.ABCMeta):
"""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/electric distribution onto a scalar phase map. :class:`~.PhaseMapper`
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:`~.FieldData` object as input and return a
:class:`~.PhaseMap` object.
"""
_log = logging.getLogger(__name__ + '.PhaseMapper')
@abc.abstractmethod
def __call__(self, field_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 field.
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 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.
"""
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:`~.VectorData`
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, magdata):
assert isinstance(magdata, VectorData), 'Only VectorData objects can be mapped!'
assert magdata.a == self.kernel.a, 'Grid spacing has to match!'
assert magdata.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert magdata.dim[1:3] == self.kernel.dim_uv, 'Dimensions do not match!'
# Process input parameters:
self.u_mag[self.kernel.slice_mag] = magdata.field[0, 0, ...] # u-component
self.v_mag[self.kernel.slice_mag] = magdata.field[1, 0, ...] # v-component
return PhaseMap(magdata.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.u_fft + self.v_mag_fft * self.kernel.v_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.u_fft
v_phase_adj_fft = mag_adj_fft * self.kernel.v_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.ravel(), -v_phase_adj.ravel()))
# TODO: Why minus?
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:`~.VectorData` 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, magdata):
self._log.debug('Calling __call__')
assert isinstance(magdata, VectorData), 'Only VectorData objects can be mapped!'
assert magdata.a == self.a, 'Grid spacing has to match!'
assert magdata.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert magdata.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
v_dim, u_dim = self.dim_uv
u_mag, v_mag = magdata.field[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.rfft2(u_mag_pad)
v_mag_fft = np.fft.rfft2(v_mag_pad)
# Calculate the Fourier transform of the phase:
f_u = np.fft.rfftfreq(u_dim + 2 * u_pad, self.a)
f_v = np.fft.fftfreq(v_dim + 2 * v_pad, self.a)
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(phase_fft)
phase = phase_pad[v_pad:v_pad + v_dim, u_pad:u_pad + u_dim]
return PhaseMap(magdata.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 = VectorData(self.a, np.zeros((3, 1) + self.dim_uv, dtype=np.float32))
magnitude_proj = np.reshape(vector, (2,) + self.dim_uv)
mag_proj.field[: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:`~.ScalarData` 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 300000.
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 = self.m
# 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 * np.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, elec_data):
self._log.debug('Calling __call__')
assert isinstance(elec_data, ScalarData), 'Only ScalarData objects can be mapped!'
assert elec_data.a == self.a, 'Grid spacing has to match!'
assert elec_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert elec_data.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
phase = self.coeff * np.squeeze(elec_data.get_mask(self.threshold))
return PhaseMap(elec_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 electrostatic field of every pixel (row-wise).
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
"""
raise NotImplementedError() # TODO: Implement right!
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 ``N**2`` entries like an electrostatic projection.
"""
raise NotImplementedError() # TODO: Implement right!
class PhaseMapperCharge(PhaseMapper):
""""""
def __init__(self, a, dim_uv, biprism_vec, v_acc=300000):
self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.biprism_vec = biprism_vec
self.v_acc = v_acc
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 * np.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 = C_e * Q_E / (4 * np.pi * EPS_0)
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, v_acc=%r)' % \
(self.__class__, self.a, self.dim_uv, self.v_acc)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMapperCharge(a=%s, dim_uv=%s, v_acc=%s)' % \
(self.a, self.dim_uv, self.v_acc)
def __call__(self, elec_data):
""" phase_dipoles() is to caculate the phase from many electric dipoles.
field include the amount of charge in every grid, unit:electron.
The normal vector of the electrode is (a,b), and the distance to the origin is the norm
of (a,b). R is the sampling rate,pixel/nm."""
R = 1 / self.a
field = elec_data.field[0, ...]
bu, bv = self.biprism_vec # biprism vector (orthogonal to biprism from origin)
bn = bu ** 2 + bv ** 2 # norm of the biprism vector
dim_v, dim_u = field.shape
# Find list of charge locations and charge values:
vq, uq = np.nonzero(field)
q = field[vq, uq]
# Calculate locations of image charges:
vm = (bu ** 2 - bv ** 2) / bn * vq - 2 * bu * bv / bn * uq + 2 * bv
um = (bv ** 2 - bu ** 2) / bn * uq - 2 * bu * bv / bn * vq + 2 * bu
# Calculate phase contribution for each charge:
phase = np.zeros((dim_v, dim_u))
for i in range(len(q)):
for u in range(dim_u):
for v in range(dim_v):
rq = np.sqrt((u - uq[i]) ** 2 + (v - vq[i]) ** 2) # charge distance
rm = np.sqrt((u - um[i]) ** 2 + (v - vm[i]) ** 2) # mirror distance
# distance
z1 = (R / 2) ** 2 - rq ** 2
z2 = (R / 2) ** 2 - rm ** 2
if z1 < 0 and z2 < 0:
phase[v, u] += -q[i] * self.coeff * np.log((rq ** 2) / (rm ** 2))
elif z1 >= 0 >= z2:
z3 = np.sqrt(z1)
phase[v, u] += (-q[i] * self.coeff *
np.log(((z3 + R / 2) ** 2) / (rm ** 2))
+ 2 * q[i] * Q_E * z3 / 2 / np.pi / EPS_0 / R)
else:
z4 = np.sqrt(z2)
phase[v, u] += (-q[i] * self.coeff *
np.log((rq ** 2) / ((z4 + R / 2) ** 2))
- 2 * q[i] * Q_E * z4 / 2 / np.pi / EPS_0 / R)
return PhaseMap(self.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 electrostatic field of every pixel (row-wise).
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
"""
raise NotImplementedError() # TODO: Implement right!
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 ``N**2`` entries like an electrostatic projection.
"""
raise NotImplementedError() # TODO: Implement right!