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 3567 deletions
pyramid/phasemap.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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 numpy as np
from PIL import Image
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.ticker import MaxNLocator
from mpl_toolkits.mplot3d import Axes3D
import cmocean
from scipy import ndimage
import warnings
from . import colors
from . import plottools
__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,
u'nrad': 1E9,
u'1/rad': 1E0,
u'1/mrad': 1E-3,
u'1/µrad': 1E-6,
u'1/nrad': 1E-9}
@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, not {}!'.format(phase.shape)
self._phase = phase.astype(dtype=np.float32)
self._dim_uv = phase.shape
@property
def phase_vec(self):
"""Vector containing the phase shift."""
return self.phase.ravel()
@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!'
confidence = confidence.astype(dtype=np.float32)
confidence /= confidence.max() + 1E-30 # Normalise!
else:
confidence = np.ones_like(self.phase, dtype=np.float32)
self._confidence = confidence
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 = np.minimum(self.confidence, other.confidence) # use minimum confidence!
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 __getitem__(self, item):
return PhaseMap(self.a, self.phase[item], self.mask[item], self.confidence[item])
def __array__(self, dtype=None): # Used for numpy ufuncs, together with __array_wrap__!
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 = ndimage.zoom(self.phase, zoom=2 ** n, order=order)
self.mask = ndimage.zoom(self.mask, zoom=2 ** n, order=0)
self.confidence = ndimage.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]]
def flip(self, axis='u'):
"""Flip/mirror the phase map around the specified axis.
Parameters
----------
axis: {'u', 'v'}, optional
The axis around which the phase map is flipped.
Returns
-------
phasemap_flip: :class:`~.PhaseMap`
A flipped copy of the :class:`~.PhaseMap` object.
"""
self._log.debug('Calling flip')
if axis == 'u':
return PhaseMap(self.a, np.flipud(self.phase), np.flipud(self.mask),
np.flipud(self.confidence))
if axis == 'v':
return PhaseMap(self.a, np.fliplr(self.phase), np.fliplr(self.mask),
np.fliplr(self.confidence))
else:
raise ValueError("Wrong input! 'u', 'v' allowed!")
def rotate(self, angle):
"""Rotate the phase map (right hand rotation).
Parameters
----------
angle: float
The angle around which the phase map is rotated.
Returns
-------
phasemap_rot: :class:`~.PhaseMap`
A rotated copy of the :class:`~.PhaseMap` object.
"""
self._log.debug('Calling rotate')
phase_rot = ndimage.rotate(self.phase, angle, reshape=False)
mask_rot = ndimage.rotate(self.mask, angle, reshape=False, order=0)
conf_rot = ndimage.rotate(self.confidence, angle, reshape=False)
return PhaseMap(self.a, phase_rot, mask_rot, conf_rot)
def shift(self, shift):
"""Shift the phase map (subpixel accuracy).
Parameters
----------
shift : float or sequence, optional
The shift along the axes. If a float, shift is the same for each axis.
If a sequence, shift should contain one value for each axis.
Returns
-------
phasemap_shift: :class:`~.PhaseMap`
A shifted copy of the :class:`~.PhaseMap` object.
"""
self._log.debug('Calling shift')
phase_rot = ndimage.shift(self.phase, shift, mode='constant', cval=0)
mask_rot = ndimage.shift(self.mask, shift, mode='constant', cval=False, order=0)
conf_rot = ndimage.shift(self.confidence, shift, mode='constant', cval=0)
return PhaseMap(self.a, phase_rot, mask_rot, conf_rot)
@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:
# noinspection PyUnresolvedReferences
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 signal:
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, unit='auto', vmin=None, vmax=None, sigma_clip=None, symmetric=True,
show_mask=True, show_conf=True, norm=None, cbar=True, # specific to plot_phase!
cmap=None, interpolation='none', axis=None, figsize=None, **kwargs):
"""Display the phasemap as a colormesh.
Parameters
----------
unit: {'rad', 'mrad', 'µrad', '1/rad', '1/mrad', '1/µrad'}, optional
The plotting unit of the phase map. The phase is scaled accordingly before plotting.
Inverse radians should be used for gain maps!
vmin : float, optional
Minimum value used for determining the plot limits. If not set, it will be
determined by the minimum of the phase directly.
vmax : float, optional
Maximum value used for determining the plot limits. If not set, it will be
determined by the minimum of the phase directly.
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`.
symmetric : boolean, optional
If True (default), a zero symmetric colormap is assumed and a zero value (which
will always be present) will be set to the central color color of the colormap.
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.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
cbar : bool, optional
If True (default), a colorbar will be plotted.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
figsize : tuple of floats (N=2)
Size of the plot figure.
Returns
-------
axis, cbar: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
Notes
-----
Uses :func:`~.plottools.format_axis` at the end. According keywords can also be given here.
"""
self._log.debug('Calling plot_phase')
a = self.a
if figsize is None:
figsize = plottools.FIGSIZE_DEFAULT
# Take units into consideration:
if unit == 'auto': # Try to automatically determine unit (recommended):
for key, value in self.UNITDICT.items():
if not key.startswith('1/'):
order = np.floor(np.log10(np.abs(self.phase).max() * value))
if -1 <= order < 2:
unit = key
if unit == 'auto': # No fitting unit was found:
unit = 'rad'
# Scale phase and make last check if order is okay:
phase = self.phase * self.UNITDICT[unit]
order = np.floor(np.log10(np.abs(phase).max()))
if order > 2 or order < -6: # Display would look bad
unit = '{} x 1E{:g}'.format(unit, order)
phase /= 10 ** order
# Calculate limits if necessary (not necessary if both limits are already set):
if vmin is None and vmax 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 limits if necessary (zero has to be present!):
if vmin is None:
vmin = np.min(phase_l)
if vmax is None:
vmax = np.max(phase_l)
# Configure colormap, to fix white to zero if colormap is symmetric:
if symmetric:
if cmap is None:
cmap = plt.get_cmap('RdBu') # TODO: use cmocean.cm.balance (flipped colours!)
# TODO: get default from "colors" or "plots" package
# TODO: make flexible, cmocean and matplotlib...
elif isinstance(cmap, str): # Get colormap if given as string:
cmap = plt.get_cmap(cmap)
vmin, vmax = np.min([vmin, -0]), np.max([0, vmax]) # Ensure zero is present!
limit = np.max(np.abs([vmin, vmax]))
start = (vmin + limit) / (2 * limit)
end = (vmax + limit) / (2 * limit)
cmap_colors = cmap(np.linspace(start, end, 256))
cmap = LinearSegmentedColormap.from_list('Symmetric', cmap_colors)
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure(figsize=figsize)
axis = fig.add_subplot(1, 1, 1)
tight = True
else:
tight = False
axis.set_aspect('equal')
# Plot the phasemap:
im = axis.imshow(phase, cmap=cmap, vmin=vmin, vmax=vmax, 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.cmaps['transparent_confidence']
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)
# Determine colorbar title:
cbar_label = kwargs.pop('cbar_label', None)
cbar_mappable = None
if cbar:
cbar_mappable = im
if cbar_label is None:
if unit.startswith('1/'):
cbar_name = 'gain'
else:
cbar_name = 'phase'
if mpl.rcParams['text.usetex'] and 'µ' in unit: # Make sure µ works in latex:
mpl.rc('text.latex', preamble=r'\usepackage{txfonts},\usepackage{lmodern}')
unit = unit.replace('µ', '$\muup$') # Upright µ!
cbar_label = u'{} [{}]'.format(cbar_name, unit)
# Return formatted axis:
return plottools.format_axis(axis, sampling=a, cbar_mappable=cbar_mappable,
cbar_label=cbar_label, tight_layout=tight, **kwargs)
def plot_holo(self, gain='auto', # specific to plot_holo!
cmap=None, interpolation='none', axis=None, figsize=None, **kwargs):
"""Display the color coded holography image.
Parameters
----------
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.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
figsize : tuple of floats (N=2)
Size of the plot figure.
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
Notes
-----
Uses :func:`~.plottools.format_axis` at the end. According keywords can also be given here.
"""
self._log.debug('Calling plot_holo')
a = self.a
if figsize is None:
figsize = plottools.FIGSIZE_DEFAULT
# Calculate gain if 'auto' is selected:
if gain == 'auto':
gain = 4 * 2 * np.pi / (np.abs(self.phase).max() + 1E-30)
gain = round(gain, -int(np.floor(np.log10(abs(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:
# 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)
# Clip outliers:
sigma_clip = 2
outlier_x = np.abs(grad_x - np.mean(grad_x)) < sigma_clip * np.std(grad_x)
grad_x_sigma = np.where(outlier_x, grad_x, np.nan)
grad_x_min, grad_x_max = np.nanmin(grad_x_sigma), np.nanmax(grad_x_sigma)
grad_x = np.clip(grad_x, grad_x_min, grad_x_max)
outlier_y = np.abs(grad_y - np.mean(grad_y)) < sigma_clip * np.std(grad_y)
grad_y_sigma = np.where(outlier_y, grad_y, np.nan)
grad_y_min, grad_y_max = np.nanmin(grad_y_sigma), np.nanmax(grad_y_sigma)
grad_y = np.clip(grad_y, grad_y_min, grad_y_max)
# Calculate colors:
if cmap is None:
cmap = colors.CMAP_CIRCULAR_DEFAULT
vector = np.asarray((grad_x, -grad_y, np.zeros_like(grad_x)))
rgb = cmap.rgb_from_vector(vector)
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(figsize=figsize)
axis = fig.add_subplot(1, 1, 1)
tight = True
else:
tight = False
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]))
note = kwargs.pop('note', None)
if note is None:
note = 'gain: {:g}'.format(gain)
stroke = kwargs.pop('stroke', 'k') # Default for holo is white with black outline!
return plottools.format_axis(axis, sampling=a, note=note, tight_layout=tight,
stroke=stroke, **kwargs)
def plot_combined(self, title='', phase_title='', holo_title='', figsize=None, **kwargs):
"""Display the phase map and the resulting color coded holography image in one plot.
Parameters
----------
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
figsize : tuple of floats (N=2)
Size of the plot figure.
Returns
-------
phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
The axes on which the graphs are plotted.
Notes
-----
Uses :func:`~.plottools.format_axis` at the end. According keywords can also be given here.
"""
self._log.debug('Calling plot_combined')
# Create combined plot and set title:
if figsize is None:
figsize = (plottools.FIGSIZE_DEFAULT[0]*2 + 1, plottools.FIGSIZE_DEFAULT[1])
fig = plt.figure(figsize=figsize)
fig.suptitle(title, fontsize=20)
# Only phase is annotated, holo will show gain:
note = kwargs.pop('note', None)
# Plot holography image:
holo_axis = fig.add_subplot(1, 2, 1, aspect='equal')
self.plot_holo(axis=holo_axis, title=holo_title, note=None, **kwargs)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
self.plot_phase(axis=phase_axis, title=phase_title, note=note, **kwargs)
# Tighten layout if axis was created here:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
plt.tight_layout()
# Return the plotting axes:
return phase_axis, holo_axis
def plot_phase_with_hist(self, bins='auto', unit='rad',
title='', phase_title='', hist_title='', figsize=None, **kwargs):
"""Display the phase map and a histogram of the phase values of all pixels.
Parameters
----------
bins : int or sequence of scalars or str, optional
Bin argument that goes to the matplotlib.hist function (more documentation there).
The default is 'auto', which tries to pick something nice.
unit: {'rad', 'mrad', 'µrad', '1/rad', '1/mrad', '1/µrad'}, optional
The plotting unit of the phase map. The phase is scaled accordingly before plotting.
Inverse radians should be used for gain maps!
title : string, optional
The super title of the plot. The default is 'Combined Plot'.
phase_title : string, optional
The title of the phase map.
hist_title : string, optional
The title of the histogram.
figsize : tuple of floats (N=2)
Size of the plot figure.
Returns
-------
phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
The axes on which the graphs are plotted.
Notes
-----
Uses :func:`~.plottools.format_axis` at the end. According keywords can also be given here.
"""
self._log.debug('Calling plot_phase_with_hist')
# Create combined plot and set title:
if figsize is None:
figsize = (plottools.FIGSIZE_DEFAULT[0]*2 + 1, plottools.FIGSIZE_DEFAULT[1])
fig = plt.figure(figsize=figsize)
fig.suptitle(title, fontsize=20)
# Plot histogram:
hist_axis = fig.add_subplot(1, 2, 1)
vec = self.phase_vec * self.UNITDICT[unit] # Take units into consideration:
vec *= np.where(self.confidence > 0.5, 1, 0).ravel() # Discard low confidence points!
hist_axis.hist(vec, bins=bins, histtype='stepfilled', color='g')
# Format histogram:
x0, x1 = hist_axis.get_xlim()
y0, y1 = hist_axis.get_ylim()
hist_axis.set(aspect=np.abs(x1 - x0) / np.abs(y1 - y0) * 0.94) # Last value because cbar!
fontsize = kwargs.get('fontsize', 16)
hist_axis.tick_params(axis='both', which='major', labelsize=fontsize)
hist_axis.set_title(hist_title, fontsize=fontsize)
hist_axis.set_xlabel('phase [{}]'.format(unit), fontsize=fontsize)
hist_axis.set_ylabel('count', fontsize=fontsize)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect=1)
self.plot_phase(unit=unit, axis=phase_axis, title=phase_title, **kwargs)
# Tighten layout:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
plt.tight_layout()
# Return the plotting axes:
return phase_axis, hist_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))
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
pyramid/phasemapper.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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 jutil 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 = np.zeros(kernel.dim_pad, dtype=kernel.u.dtype)
self.v_mag = np.zeros(kernel.dim_pad, dtype=kernel.u.dtype)
self.phase_adj = np.zeros(kernel.dim_pad, dtype=kernel.u.dtype)
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.phase_adj[self.kernel.slice_phase] = vector.reshape(self.kernel.dim_uv)
phase_adj_fft = fft.irfft2_adj(self.phase_adj)
u_mag_adj_fft = phase_adj_fft * np.conj(self.kernel.u_fft)
v_mag_adj_fft = phase_adj_fft * np.conj(self.kernel.v_fft)
u_mag_adj = fft.rfft2_adj(u_mag_adj_fft)[self.kernel.slice_mag]
v_mag_adj = fft.rfft2_adj(v_mag_adj_fft)[self.kernel.slice_mag]
result = np.concatenate((u_mag_adj.ravel(), v_mag_adj.ravel()))
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):
"""""" # TODO: Write Docstring!
def __init__(self, a, dim_uv, electrode_vec, v_acc=300000):
self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.electrode_vec = electrode_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):
""" This is to calculate the phase from many electric dipoles. The model used to avoid singularity is
the homogeneously-distributed charged metallic sphere.
The elec_data includes the amount of charge in every grid, unit:electron.
The electrode_vec is the normal vector of the electrode, (elec_a,elec_b), and the distance to the origin is
the norm of (elec_a,elec_b).
R_sam is the sampling rate,pixel/nm."""
R_sam = 1 / self.a
R = R_sam / 2 # The radius of the charged sphere
field = elec_data.field[0, ...]
elec_a, elec_b = self.electrode_vec # electrode vector (orthogonal to the electrode from origin)
elec_n = elec_a ** 2 + elec_b ** 2
dim_v, dim_u = field.shape
u_cor = range(0, dim_u, 1)
v_cor = range(0, dim_v, 1)
v, u = np.meshgrid(v_cor, u_cor)
# 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
um = (elec_b ** 2 - elec_a ** 2) / elec_n * uq - 2 * elec_a * elec_b / elec_n * vq + 2 * elec_a
vm = (elec_a ** 2 - elec_b ** 2) / elec_n * vq - 2 * elec_a * elec_b / elec_n * uq + 2 * elec_b
# Calculate phase contribution for each charge:
phase = np.zeros((dim_v, dim_u))
for i in range(len(q)):
# The projected distance from the charges or image charges
r1 = np.sqrt((u - uq[i]) ** 2 + (v - vq[i]) ** 2)
r2 = np.sqrt((u - um[i]) ** 2 + (v - vm[i]) ** 2)
# The square height when the path come across the sphere
z1 = R ** 2 - r1 ** 2
z2 = R ** 2 - r2 ** 2
# Phase calculation in 3 different cases
# case 1 totally out of the sphere
case1 = ((z1 < 0) & (z2 < 0))
phase[case1] += - q[i] * self.coeff * np.log((r1[case1] ** 2) / (r2[case1] ** 2))
# case 2: inside the charge sphere
case2 = ((z1 >= 0) & (z2 <= 0))
# The height when the path come across the charge sphere
z3 = np.sqrt(z1)
phase[case2] += q[i] * self.coeff * (- np.log((z3[case2] + R) ** 2 / r2[case2] ** 2) +
(2 * z3[case2] / R + 2 * z3[case2] ** 3 / 3 / R ** 3))
# case 3 : inside the image charge sphere
case3 = np.logical_not(case1 | case2)
# The height whe the path comes across the image charge sphere
z4 = np.sqrt(z2)
phase[case3] += q[i] * self.coeff * (np.log((z4[case3] + R) ** 2 / r1[case3] ** 2) -
(2 * z4[case3] / R + 2 * z4[case3] ** 3 / 3 / R ** 3))
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!
pyramid/plottools.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
# Adapted from mpl_toolkits.axes_grid2
"""This module provides the useful plotting utilities."""
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.offsetbox import AnchoredOffsetbox, AuxTransformBox, VPacker, TextArea
from matplotlib.transforms import blended_transform_factory, IdentityTransform
from matplotlib.patches import Rectangle
from matplotlib import patheffects
from matplotlib.ticker import MaxNLocator, FuncFormatter
from mpl_toolkits.axes_grid.inset_locator import inset_axes
from mpl_toolkits.axes_grid1 import make_axes_locatable
import warnings
from . import colors
__all__ = ['format_axis', 'pretty_plots', 'add_scalebar',
'add_annotation', 'add_colorwheel', 'add_cbar']
FIGSIZE_DEFAULT = (6.7, 5) # TODO: Apparently does not fit as well as before...
FONTSIZE_DEFAULT = 20
STROKE_DEFAULT = None
def pretty_plots(figsize=None, fontsize=None, stroke=None):
"""Set IPython formats (for interactive and PDF output) and set pretty matplotlib font."""
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png', 'pdf') # png for interactive, pdf, for PDF output!
mpl.rcParams['mathtext.fontset'] = 'stix' # Mathtext in $...$!
mpl.rcParams['font.family'] = 'STIXGeneral' # Set normal text to look the same!
mpl.rcParams['figure.max_open_warning'] = 0 # Disable Max Open Figure warning!
if figsize is not None:
global FIGSIZE_DEFAULT
FIGSIZE_DEFAULT = figsize
mpl.rcParams['figure.figsize'] = FIGSIZE_DEFAULT
if fontsize is not None:
global FONTSIZE_DEFAULT
FONTSIZE_DEFAULT = fontsize
mpl.rcParams['xtick.labelsize'] = FONTSIZE_DEFAULT
mpl.rcParams['ytick.labelsize'] = FONTSIZE_DEFAULT
mpl.rcParams['axes.labelsize'] = FONTSIZE_DEFAULT
mpl.rcParams['legend.fontsize'] = FONTSIZE_DEFAULT
global STROKE_DEFAULT
STROKE_DEFAULT = stroke
def add_scalebar(axis, sampling=1, fontsize=None, stroke=None):
"""Add a scalebar to the axis.
Parameters
----------
axis : :class:`~matplotlib.axes.AxesSubplot`
Axis to which the scalebar is added.
sampling : float, optional
The grid spacing in nm. If not given, 1 nm is assumed.
fontsize : int, optional
The fontsize which should be used for the label. Default is 16.
stroke : None or color, optional
If not None, a stroke will be applied to the text, e.g. to make it more visible.
Returns
-------
aoffbox : :class:`~matplotlib.offsetbox.AnchoredOffsetbox`
The box containing the scalebar.
"""
# TODO: Background (black outline) not visible...
if fontsize is None:
fontsize = FONTSIZE_DEFAULT
if stroke is None:
stroke = STROKE_DEFAULT
# Transform that scales the width along the data and leaves height constant at 8 pt (text):
transform = blended_transform_factory(axis.transData, IdentityTransform())
# Transform axis borders (1, 1) to data borders to get number of pixels in y and x:
bb0 = axis.transLimits.inverted().transform((0, 0))
bb1 = axis.transLimits.inverted().transform((1, 1))
dim_uv = (int(abs(bb1[1] - bb0[1])), int(abs(bb1[0] - bb0[0])))
# Calculate scale
scale = np.max((dim_uv[1] / 4, 1)) * sampling
thresholds = [1, 5, 10, 50, 100, 500, 1000]
for t in thresholds: # For larger grids (real images), multiples of threshold look better!
if scale > t:
scale = (scale // t) * t
# Set dimensions:
width = scale / sampling # In data coordinate system!
height = 8 # In display coordinate system!
# Set label:
if scale >= 1000: # Use higher order instead!
label = '{:.3g} µm'.format(scale/1000)
else:
label = '{:.3g} nm'.format(scale)
# Create scalebar rectangle:
bars = AuxTransformBox(transform)
bars.add_artist(Rectangle((0, 0), width, height, fc='w', linewidth=1,
clip_box=axis.bbox, clip_on=True))
# Create text:
txtcolor = 'w' if stroke == 'k' else 'k'
txt = TextArea(label, textprops={'color': txtcolor, 'fontsize': fontsize})
txt.set_clip_box(axis.bbox)
if stroke is not None:
txt._text.set_path_effects([patheffects.withStroke(linewidth=2, foreground=stroke)])
# Pack both together an create AnchoredOffsetBox:
bars = VPacker(children=[txt, bars], align="center", pad=0.5, sep=5)
aoffbox = AnchoredOffsetbox(loc=3, pad=0.5, borderpad=0.1, child=bars, frameon=False)
axis.add_artist(aoffbox)
# Return:
return aoffbox
def add_annotation(axis, label, stroke=None, fontsize=None):
"""Add an annotation to the axis on the upper left corner.
Parameters
----------
axis : :class:`~matplotlib.axes.AxesSubplot`
Axis to which the annotation is added.
label : string
The text of the annotation.
fontsize : int, optional
The fontsize which should be used for the annotation. Default is 16.
stroke : None or color, optional
If not None, a stroke will be applied to the text, e.g. to make it more visible.
Returns
-------
aoffbox : :class:`~matplotlib.offsetbox.AnchoredOffsetbox`
The box containing the annotation.
"""
if fontsize is None:
fontsize = FONTSIZE_DEFAULT
if stroke is None:
stroke = STROKE_DEFAULT
# Create text:
txtcolor = 'w' if stroke == 'k' else 'k'
txt = TextArea(label, textprops={'color': txtcolor, 'fontsize': fontsize})
txt.set_clip_box(axis.bbox)
if stroke is not None:
txt._text.set_path_effects([patheffects.withStroke(linewidth=2, foreground=stroke)])
# Pack into and add AnchoredOffsetBox:
aoffbox = AnchoredOffsetbox(loc=2, pad=0.5, borderpad=0.1, child=txt, frameon=False)
axis.add_artist(aoffbox)
return aoffbox
def add_colorwheel(axis):
"""Add a colorwheel to the axis on the upper right corner.
Parameters
----------
axis : :class:`~matplotlib.axes.AxesSubplot`
Axis to which the colorwheel is added.
Returns
-------
axis : :class:`~matplotlib.axes.AxesSubplot`
The same axis which was given to this function is returned.
"""
from mpl_toolkits.axes_grid.inset_locator import inset_axes
inset_axes = inset_axes(axis, width=0.75, height=0.75, loc=1)
inset_axes.axis('off')
cmap = colors.CMAP_CIRCULAR_DEFAULT
bgcolor = axis.get_facecolor()
return cmap.plot_colorwheel(size=100, axis=inset_axes, alpha=0, bgcolor=bgcolor, arrows=True)
def add_cbar(axis, mappable, label='', fontsize=None):
"""Add a colorbar to the right of the given axis.
Parameters
----------
axis : :class:`~matplotlib.axes.AxesSubplot`
Axis to which the colorbar is added.
mappable : mappable pyplot
If this is not None, a colorbar will be plotted on the right side of the axes,
label : string, optional
The label of the colorbar. If not set, no label is used.
fontsize : int, optional
The fontsize which should be used for the label. Default is 16.
Returns
-------
cbar : :class:`~matplotlib.Colorbar`
The created colorbar.
"""
if fontsize is None:
fontsize = FONTSIZE_DEFAULT
divider = make_axes_locatable(axis)
cbar_ax = divider.append_axes('right', size='5%', pad=0.1)
cbar = plt.colorbar(mappable, cax=cbar_ax)
cbar.ax.tick_params(labelsize=fontsize)
# Make sure labels don't stick out of tight bbox:
labels = cbar.ax.get_yticklabels()
delta = 0.03 * (cbar.vmax - cbar.vmin)
lmin = float(labels[0]._text.replace(u'\u2212', '-').strip('$')) # No unicode or latex!
lmax = float(labels[-1]._text.replace(u'\u2212', '-').strip('$')) # No unicode or latex!
redo_max = True if cbar.vmax - lmax < delta else False
redo_min = True if lmin - cbar.vmin < delta else False
mappable.set_clim(cbar.vmin - delta * redo_min, cbar.vmax + delta * redo_max)
# A lot of plotting magic to make plot width consistent (mostly):
cbar.ax.set_yticklabels(labels, ha='right')
renderer = plt.gcf().canvas.get_renderer()
bbox_0 = labels[0].get_window_extent(renderer)
bbox_1 = labels[-1].get_window_extent(renderer)
cbar_pad = np.max((bbox_0.width, bbox_1.width)) + 5 # bit of padding left!
cbar.ax.yaxis.set_tick_params(pad=cbar_pad)
max_txt = plt.text(0, 0, u'\u22120.00', fontsize=fontsize)
bbox_max = max_txt.get_window_extent(renderer)
max_txt.remove()
cbar_pad_max = bbox_max.width + 10 # bit of padding right!
cbar.set_label(label, fontsize=fontsize, labelpad=max(cbar_pad_max - cbar_pad, 0))
# Set focus back to axis and return cbar:
plt.sca(axis)
return cbar
def add_coords(axis, coords=('x', 'y')):
ins_ax = inset_axes(axis, width="5%", height="5%", loc=3, borderpad=2.2)
if coords == 3:
coords = ('x', 'y', 'z')
elif coords == 2:
coords = ('x', 'y')
if len(coords) == 3:
ins_ax.arrow(0.5, 0.45, -1.05, -0.75, fc="k", ec="k",
head_width=0.2, head_length=0.3, linewidth=3, clip_on=False)
ins_ax.arrow(0.5, 0.45, 0.96, -0.75, fc="k", ec="k",
head_width=0.2, head_length=0.3, linewidth=3, clip_on=False)
ins_ax.arrow(0.5, 0.45, 0, 1.35, fc="k", ec="k",
head_width=0.2, head_length=0.3, linewidth=3, clip_on=False)
ins_ax.annotate(coords[0], xy=(0, 0), xytext=(-0.9, 0), fontsize=20, clip_on=False)
ins_ax.annotate(coords[1], xy=(0, 0), xytext=(1.4, 0.1), fontsize=20, clip_on=False)
ins_ax.annotate(coords[2], xy=(0, 0), xytext=(0.7, 1.5), fontsize=20, clip_on=False)
elif len(coords) == 2:
ins_ax.arrow(-0.5, -0.5, 1.5, 0, fc="k", ec="k",
head_width=0.2, head_length=0.3, linewidth=3, clip_on=False)
ins_ax.arrow(-0.5, -0.5, 0, 1.5, fc="k", ec="k",
head_width=0.2, head_length=0.3, linewidth=3, clip_on=False)
ins_ax.annotate(coords[0], xy=(0, 0), xytext=(1.3, -0.05), fontsize=20, clip_on=False)
ins_ax.annotate(coords[1], xy=(0, 0), xytext=(-0.2, 1.3), fontsize=20, clip_on=False)
ins_ax.axis('off')
plt.sca(axis)
# TODO: These parameters in other plot functions belong in a dedicated dictionary!!!
def format_axis(axis, format_axis=True, title='', fontsize=None, stroke=None, scalebar=True,
hideaxes=None, sampling=1, note=None, colorwheel=False, cbar_mappable=None,
cbar_label='', tight_layout=True, keep_labels=False, coords=None, **_):
"""Format an axis and add a lot of nice features.
Parameters
----------
axis : :class:`~matplotlib.axes.AxesSubplot`
Axis on which the graph is plotted.
format_axis : bool, optional
If False, the formatting will be skipped (the axis is still returned). Default is True.
title : string, optional
The title of the plot. The default is an empty string''.
fontsize : int, optional
The fontsize which should be used for labels and titles. Default is 16.
stroke : None or color, optional
If not None, a stroke will be applied to the text, e.g. to make it more visible.
scalebar : bool, optional
Defines if a scalebar should be plotted in the lower left corner (default: True). Axes
are made invisible. If set to False, the axes are formatted to ook pretty, instead.
hideaxes : True, optional
If True, the axes will be turned invisible. If not specified (None), this is True if a
scalebar is plotted, False otherwise.
sampling : float, optional
The grid spacing in nm. If not given, 1 nm is assumed.
note: string or None, optional
An annotation string which is displayed in the upper left
colorwheel : bool, optional
Defines if a colorwheel should be plotted in the upper right corner (default: False).
cbar_mappable : mappable pyplot or None, optional
If this is not None, a colorbar will be plotted on the right side of the axes,
which uses this mappable object.
cbar_label : string, optional
The label of the colorbar. If `None`, no label is used.
tight_layout : bool, optional
If True, `plt.tight_layout()` is executed after the formatting. Default is True.
Returns
-------
axis : :class:`~matplotlib.axes.AxesSubplot`
The same axis which was given to this function is returned.
"""
if not format_axis: # Skip (sometimes useful if more than one plot is used on the same axis!
return axis
if fontsize is None:
fontsize = FONTSIZE_DEFAULT
if stroke is None:
stroke = STROKE_DEFAULT
# Add scalebar:
if scalebar:
add_scalebar(axis, sampling=sampling, fontsize=fontsize, stroke=stroke)
if hideaxes is None:
hideaxes = True
if not keep_labels:
# Set_title
axis.set_title(title, fontsize=fontsize)
# Get dimensions:
bb0 = axis.transLimits.inverted().transform((0, 0))
bb1 = axis.transLimits.inverted().transform((1, 1))
dim_uv = (int(abs(bb1[1] - bb0[1])), int(abs(bb1[0] - bb0[0])))
# Set the title and the axes labels:
axis.set_xlim(0, dim_uv[1])
axis.set_ylim(0, dim_uv[0])
# Determine major tick locations (useful for grid, even if ticks will not be used):
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))
# Hide axes label and ticks if wanted:
if hideaxes:
for tic in axis.xaxis.get_major_ticks():
tic.tick1On = tic.tick2On = False
tic.label1On = tic.label2On = False
for tic in axis.yaxis.get_major_ticks():
tic.tick1On = tic.tick2On = False
tic.label1On = tic.label2On = False
else: # Set the axes ticks and labels:
if not keep_labels:
axis.set_xlabel('u-axis [nm]')
axis.set_ylabel('v-axis [nm]')
func_formatter = FuncFormatter(lambda x, pos: '{:.3g}'.format(x * sampling))
axis.xaxis.set_major_formatter(func_formatter)
axis.yaxis.set_major_formatter(func_formatter)
axis.tick_params(axis='both', which='major', labelsize=fontsize)
axis.xaxis.label.set_size(fontsize)
axis.yaxis.label.set_size(fontsize)
# Add annotation:
if note:
add_annotation(axis, label=note, fontsize=fontsize, stroke=stroke)
# Add coords:
if coords:
add_coords(axis, coords=coords)
# Add colorhweel:
if colorwheel:
add_colorwheel(axis)
# Add colorbar:
if cbar_mappable:
# Construct colorbar:
add_cbar(axis, mappable=cbar_mappable, label=cbar_label, fontsize=fontsize)
# Tighten layout if axis was created here:
if tight_layout:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
plt.tight_layout()
# Return plotting axis:
return axis
pyramid/projector.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the abstract base class :class:`~.Projector` and concrete subclasses for
projections of vector and scalar fields."""
import itertools
import logging
try:
if type(get_ipython()).__name__ == 'ZMQInteractiveShell': # IPython Notebook!
from tqdm import tqdm_notebook as tqdm
else: # IPython, but not a Notebook (e.g. terminal)
from tqdm import tqdm
except NameError:
from tqdm import tqdm
import numpy as np
from numpy import pi
from scipy.sparse import coo_matrix, csr_matrix
from pyramid.fielddata import VectorData, ScalarData
from pyramid.quaternion import Quaternion
__all__ = ['RotTiltProjector', 'XTiltProjector', 'YTiltProjector', 'SimpleProjector']
class Projector(object):
"""Base class representing a projection function.
The :class:`~.Projector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid. :class:`~.Projector` is an abstract base
class and provides a unified interface which should be subclassed with a custom
:func:`__init__` function, which should call the parent :func:`__init__` method. Concrete
subclasses can be called as a function and take a `vector` as argument which contains the
3-dimensional field. The output is the projected field, given as a `vector`. Depending on the
length of the input and the given dimensions `dim` at construction time, vector or scalar
projection is choosen intelligently.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
dim_uv : tuple (N=2)
Dimensions (v, u) of the projected grid.
size_3d : int
Number of voxels of the 3-dimensional grid.
size_2d : int
Number of pixels of the 2-dimensional projected grid.
weight : :class:`~scipy.sparse.csr_matrix` (N=2)
The weight matrix containing the weighting coefficients for the 3D to 2D mapping.
coeff : list (N=2)
List containing the six weighting coefficients describing the influence of the 3 components
of a 3-dimensional vector field on the 2 projected components.
m: int
Size of the image space.
n: int
Size of the input space.
sparsity : float
Measures the sparsity of the weighting (not the complete one!), 1 means completely sparse!
"""
_log = logging.getLogger(__name__ + '.Projector')
@property
def sparsity(self):
"""The sparsity of the projector weight matrix."""
return 1. - len(self.weight.data) / np.prod(self.weight.shape)
def __init__(self, dim, dim_uv, weight, coeff):
self._log.debug('Calling __init__')
self.dim = tuple(dim)
self.dim_uv = tuple(dim_uv)
self.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
self.n = 3 * np.prod(dim)
self.m = 2 * np.prod(dim_uv)
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(dim=%r, dim_uv=%r, weight=%r, coeff=%r)' % \
(self.__class__, self.dim, self.dim_uv, self.weight, self.coeff)
def __str__(self):
self._log.debug('Calling __str__')
return 'Projector(dim=%s, dim_uv=%s, coeff=%s)' % (self.dim, self.dim_uv, self.coeff)
def __call__(self, field_data):
if isinstance(field_data, VectorData):
field_empty = np.zeros((3, 1) + self.dim_uv, dtype=field_data.field.dtype)
field_data_proj = VectorData(field_data.a, field_empty)
field_proj = self.jac_dot(field_data.field_vec).reshape((2,) + self.dim_uv)
field_data_proj.field[0:2, 0, ...] = field_proj
elif isinstance(field_data, ScalarData):
field_empty = np.zeros((1,) + self.dim_uv, dtype=field_data.field.dtype)
field_data_proj = ScalarData(field_data.a, field_empty)
field_proj = self.jac_dot(field_data.field_vec).reshape(self.dim_uv)
field_data_proj.field[0, ...] = field_proj
else:
raise TypeError('Input is neither of type VectorData or ScalarData')
return field_data_proj
def _vector_field_projection(self, vector):
result = np.zeros(2 * self.size_2d, dtype=vector.dtype)
# Go over all possible component projections (z, y, x) to (u, v):
vec_x, vec_y, vec_z = np.split(vector, 3)
vec_x_weighted = self.weight.dot(vec_x)
vec_y_weighted = self.weight.dot(vec_y)
vec_z_weighted = self.weight.dot(vec_z)
slice_u = slice(0, self.size_2d)
slice_v = slice(self.size_2d, 2 * self.size_2d)
if self.coeff[0][0] != 0: # x to u
result[slice_u] += self.coeff[0][0] * vec_x_weighted
if self.coeff[0][1] != 0: # y to u
result[slice_u] += self.coeff[0][1] * vec_y_weighted
if self.coeff[0][2] != 0: # z to u
result[slice_u] += self.coeff[0][2] * vec_z_weighted
if self.coeff[1][0] != 0: # x to v
result[slice_v] += self.coeff[1][0] * vec_x_weighted
if self.coeff[1][1] != 0: # y to v
result[slice_v] += self.coeff[1][1] * vec_y_weighted
if self.coeff[1][2] != 0: # z to v
result[slice_v] += self.coeff[1][2] * vec_z_weighted
return result
def _vector_field_projection_T(self, vector):
result = np.zeros(3 * self.size_3d)
# Go over all possible component projections (u, v) to (z, y, x):
vec_u, vec_v = np.split(vector, 2)
vec_u_weighted = self.weight.T.dot(vec_u)
vec_v_weighted = self.weight.T.dot(vec_v)
slice_x = slice(0, self.size_3d)
slice_y = slice(self.size_3d, 2 * self.size_3d)
slice_z = slice(2 * self.size_3d, 3 * self.size_3d)
if self.coeff[0][0] != 0: # u to x
result[slice_x] += self.coeff[0][0] * vec_u_weighted
if self.coeff[0][1] != 0: # u to y
result[slice_y] += self.coeff[0][1] * vec_u_weighted
if self.coeff[0][2] != 0: # u to z
result[slice_z] += self.coeff[0][2] * vec_u_weighted
if self.coeff[1][0] != 0: # v to x
result[slice_x] += self.coeff[1][0] * vec_v_weighted
if self.coeff[1][1] != 0: # v to y
result[slice_y] += self.coeff[1][1] * vec_v_weighted
if self.coeff[1][2] != 0: # v to z
result[slice_z] += self.coeff[1][2] * vec_v_weighted
return result
def _scalar_field_projection(self, vector):
self._log.debug('Calling _scalar_field_projection')
return np.array(self.weight.dot(vector))
def _scalar_field_projection_T(self, vector):
self._log.debug('Calling _scalar_field_projection_T')
return np.array(self.weight.T.dot(vector))
def jac_dot(self, vector):
"""Multiply a `vector` with the jacobi matrix of this :class:`~.Projector` object.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector containing the field which should be projected. Must have the same or 3 times
the size of `size_3d` of the projector for scalar and vector projection, respectively.
Returns
-------
proj_vector : :class:`~numpy.ndarray` (N=1)
Vector containing the projected field of the 2-dimensional grid. The length is
always`size_2d`.
"""
if len(vector) == 3 * self.size_3d: # mode == 'vector'
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for '
'vector- or scalar-field-projection!')
def jac_T_dot(self, vector):
"""Multiply a `vector` with the transp. jacobi matrix of this :class:`~.Projector` object.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector containing the field which should be projected. Must have the same or 2 times
the size of `size_2d` of the projector for scalar and vector projection, respectively.
Returns
-------
proj_vector : :class:`~numpy.ndarray` (N=1)
Vector containing the multiplication of the input with the transposed jacobi matrix
of the :class:`~.Projector` object.
"""
if len(vector) == 2 * self.size_2d: # mode == 'vector'
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for '
'vector- or scalar-field-projection!')
def save(self, filename, overwrite=True):
"""Saves the projector as an HDF5 file.
Parameters
----------
filename: str
Name of the file which the phasemap is saved into. HDF5 files are supported.
overwrite: bool, optional
If True (default), an existing file will be overwritten, if False, this
(silently!) does nothing.
"""
from .file_io.io_projector import save_projector
save_projector(self, filename, overwrite)
def get_info(self, verbose):
"""Get specific information about the projector as a string.
Parameters
----------
verbose: boolean, optional
If this is true, the text looks prettier (maybe using latex). Default is False for the
use in file names and such.
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
"""
return 'Base projector'
class RotTiltProjector(Projector):
"""Class representing a projection function with a rotation around z followed by tilt around x.
The :class:`~.XTiltProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
rotation : float
Angle in `rad` describing the rotation around the z-axis before the tilt is happening.
tilt : float
Angle in `rad` describing the tilt of the beam direction relative to the x-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set defaults to the (y, x)-dimensions.
subcount : int (optional)
Number of subpixels along one axis. This is used to create the lookup table which uses
a discrete subgrid to estimate the impact point of a voxel onto a pixel and the weight on
all surrounding pixels. Default is 11 (odd numbers provide a symmetric center).
"""
_log = logging.getLogger(__name__ + '.RotTiltProjector')
def __init__(self, dim, rotation, tilt, dim_uv=None, subcount=11, verbose=False):
self._log.debug('Calling __init__')
self.rotation = rotation
self.tilt = tilt
# Determine dimensions:
dim_z, dim_y, dim_x = dim
center = (dim_z / 2., dim_y / 2., dim_x / 2.)
if dim_uv is None:
dim_v = max(dim_x, dim_y) # first rotate around z-axis (take x and y into account)
dim_u = max(dim_v, dim_z) # then tilt around x-axis (now z matters, too)
dim_uv = (dim_v, dim_u)
dim_v, dim_u = dim_uv
# Creating coordinate list of all voxels:
voxels = list(itertools.product(range(dim_z), range(dim_y), range(dim_x)))
# Calculate vectors to voxels relative to rotation center:
voxel_vecs = (np.asarray(voxels) + 0.5 - np.asarray(center)).T
# Create tilt, rotation and combined quaternion, careful: Quaternion(w,x,y,z), not (z,y,x):
quat_x = Quaternion.from_axisangle((1, 0, 0), tilt) # Tilt around x-axis
quat_z = Quaternion.from_axisangle((0, 0, 1), rotation) # Rotate around z-axis
quat = quat_x * quat_z # Combined quaternion (first rotate around z, then tilt around x)
# Calculate impact positions on the projected pixel coordinate grid (flip because quat.):
impacts = np.flipud(quat.matrix[:2, :].dot(np.flipud(voxel_vecs))) # only care for x/y
impacts[1, :] += dim_u / 2. # Shift back to normal indices
impacts[0, :] += dim_v / 2. # Shift back to normal indices
# Calculate equivalence radius:
R = (3 / (4 * np.pi)) ** (1 / 3.)
# Prepare weight matrix calculation:
rows = [] # 2D projection
columns = [] # 3D distribution
data = [] # weights
# Create 4D lookup table (1&2: which neighbour weight?, 3&4: which subpixel is hit?)
weight_lookup = self._create_weight_lookup(subcount, R)
# Go over all voxels:
disable = not verbose
for i, voxel in enumerate(tqdm(voxels, disable=disable, leave=False,
desc='Set up projector')):
column_index = voxel[0] * dim_y * dim_x + voxel[1] * dim_x + voxel[2]
remainder, impact = np.modf(impacts[:, i]) # split index of impact and remainder!
sub_pixel = (remainder * subcount).astype(dtype=np.int) # sub_pixel inside impact px.
# Go over all influenced pixels (impact and neighbours, indices are [0, 1, 2]!):
for px_ind in list(itertools.product(range(3), range(3))):
# Pixel indices influenced by the impact (px_ind-1 to center them around impact):
pixel = (impact + np.array(px_ind) - 1).astype(dtype=np.int)
# Check if pixel is out of bound:
if 0 <= pixel[0] < dim_uv[0] and 0 <= pixel[1] < dim_uv[1]:
# Lookup weight in 4-dimensional lookup table!
weight = weight_lookup[px_ind[0], px_ind[1], sub_pixel[0], sub_pixel[1]]
# Only write into sparse matrix if weight is not zero:
if weight != 0.:
row_index = pixel[0] * dim_u + pixel[1]
columns.append(column_index)
rows.append(row_index)
data.append(weight)
# Calculate weight matrix and coefficients for jacobi matrix:
shape = (np.prod(dim_uv), np.prod(dim))
weights = csr_matrix(coo_matrix((data, (rows, columns)), shape=shape))
# Calculate coefficients by rotating unity matrix (unit vectors, (x,y,z)):
coeff = quat.matrix[:2, :].dot(np.eye(3))
super().__init__(dim, dim_uv, weights, coeff)
self._log.debug('Created ' + str(self))
@staticmethod
def _create_weight_lookup(subcount, R):
s = subcount
Rz = R * s # Radius in subgrid units
dim_zoom = (3 * s, 3 * s) # Dimensions of the subgrid, (3, 3) because of neighbour count!
cent_zoom = (np.asarray(dim_zoom) / 2.).astype(dtype=np.int) # Center of the subgrid
y, x = np.indices(dim_zoom)
y -= cent_zoom[0]
x -= cent_zoom[1]
# Calculate projected thickness of an equivalence sphere (normed!):
d = np.where(np.hypot(x, y) <= Rz, Rz ** 2 - x ** 2 - y ** 2, 0)
d = np.sqrt(d)
d /= d.sum()
# Create lookup table (4D):
lookup = np.zeros((3, 3, s, s))
# Go over all 9 pixels (center and neighbours):
for pixel in list(itertools.product(range(3), range(3))):
pixel_lb = np.array(pixel) * s # Convert to subgrid, hit bottom left of the pixel!
# Go over all subpixels in the center that can be hit:
for sub_pixel in list(itertools.product(range(s), range(s))):
shift = np.array(sub_pixel) - np.array((s // 2, s // 2)) # relative to center!
lb = pixel_lb - shift # Shift summing zone according to hit subpixel!
# Make sure, that the summing zone is in bounds (otherwise correct accordingly):
lb = np.where(lb >= 0, lb, [0, 0])
tr = np.where(lb < 3 * s, lb + np.array((s, s)), [3 * s, 3 * s])
# Calculate weight by summing over the summing zone:
weight = d[lb[0]:tr[0], lb[1]:tr[1]].sum()
lookup[pixel[0], pixel[1], sub_pixel[0], sub_pixel[1]] = weight
return lookup
def get_info(self, verbose=False):
"""Get specific information about the projector as a string.
Parameters
----------
verbose: boolean, optional
If this is true, the text looks prettier (maybe using latex). Default is False for the
use in file names and such.
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
"""
theta_ang = int(np.round(self.rotation * 180 / pi))
phi_ang = int(np.round(self.tilt * 180 / pi))
if verbose:
return u'$\\theta = {:d}$°, $\phi = {:d}$°'.format(theta_ang, phi_ang)
else:
return u'theta={:d}_phi={:d}°'.format(theta_ang, phi_ang)
class XTiltProjector(Projector):
"""Class representing a projection function with a tilt around the x-axis.
The :class:`~.XTiltProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
tilt : float
Angle in `rad` describing the tilt of the beam direction relative to the x-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set defaults to the (y, x)-dimensions.
"""
_log = logging.getLogger(__name__ + '.XTiltProjector')
def __init__(self, dim, tilt, dim_uv=None, verbose=False):
self._log.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), perpendicular (perp, y) and rotation (rot, x) axis:
dim_proj, dim_perp, dim_rot = dim
if dim_uv is None:
dim_uv = (max(dim_perp, dim_proj), dim_rot) # x-y-plane
dim_v, dim_u = dim_uv # y, x
assert dim_v >= dim_perp and dim_u >= dim_rot, 'Projected dimensions are too small!'
# Creating coordinate list of all voxels (for one slice):
voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-y-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj / 2., dim_perp / 2.)
positions = self._get_position(voxels, center, tilt, dim_v)
# Calculate weight-matrix:
r = 1 / np.sqrt(np.pi) # radius of the voxel circle
rho = 0.5 / r
row = []
col = []
data = []
# One slice:
disable = not verbose
for i, voxel in enumerate(tqdm(voxels, disable=disable, leave=False,
desc='Set up projector')):
impacts = self._get_impact(positions[i], r, dim_v) # impact along projected y-axis
voxel_index = voxel[0] * dim_rot * dim_perp + voxel[1] * dim_rot # 0: z, 1: y
for impact in impacts:
impact_index = impact * dim_u + (dim_u - dim_rot) // 2
distance = np.abs(impact + 0.5 - positions[i])
delta = distance / r
col.append(voxel_index)
row.append(impact_index)
data.append(self._get_weight(delta, rho))
# All other slices (along x):
data = np.tile(data, dim_rot)
columns = np.tile(col, dim_rot)
rows = np.tile(row, dim_rot)
addition = np.repeat(np.arange(dim_rot), len(row))
columns += addition
rows += addition
# Calculate weight matrix and coefficients for jacobi matrix:
shape = (np.prod(dim_uv), np.prod(dim))
weight = csr_matrix(coo_matrix((data, (rows, columns)), shape=shape))
coeff = [[1, 0, 0], [0, np.cos(tilt), np.sin(tilt)]]
super().__init__(dim, dim_uv, weight, coeff)
self._log.debug('Created ' + str(self))
@staticmethod
def _get_position(points, center, tilt, size):
point_vecs = np.asarray(points) + 0.5 - np.asarray(center) # vectors pointing to points
direc_vec = np.array((np.cos(tilt), -np.sin(tilt))) # vector pointing along projection
distances = np.cross(direc_vec, point_vecs) # here (special case): divisor is one!
distances += size / 2. # Shift to the center of the projection
return distances
@staticmethod
def _get_impact(pos, r, size):
return [x for x in np.arange(np.floor(pos - r), np.floor(pos + r) + 1, dtype=int)
if 0 <= x < size]
@staticmethod
def _get_weight(delta, rho): # use circles to represent the voxels
lo, up = delta - rho, delta + rho
# Upper boundary:
if up >= 1:
w_up = 0.5
else:
w_up = (up * np.sqrt(1 - up ** 2) + np.arctan(up / np.sqrt(1 - up ** 2))) / pi
# Lower boundary:
if lo <= -1:
w_lo = -0.5
else:
w_lo = (lo * np.sqrt(1 - lo ** 2) + np.arctan(lo / np.sqrt(1 - lo ** 2))) / pi
return w_up - w_lo
def get_info(self, verbose=False):
"""Get specific information about the projector as a string.
Parameters
----------
verbose: boolean, optional
If this is true, the text looks prettier (maybe using latex). Default is False for the
use in file names and such.
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
"""
if verbose:
return u'x-tilt: $\phi = {:d}$°'.format(int(np.round(self.tilt * 180 / pi)))
else:
return u'xtilt_phi={:d}°'.format(int(np.round(self.tilt * 180 / pi)))
class YTiltProjector(Projector):
"""Class representing a projection function with a tilt around the y-axis.
The :class:`~.YTiltProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
tilt : float
Angle in `rad` describing the tilt of the beam direction relative to the y-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set defaults to the (y, x)-dimensions.
"""
_log = logging.getLogger(__name__ + '.YTiltProjector')
def __init__(self, dim, tilt, dim_uv=None, verbose=False):
self._log.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
dim_proj, dim_rot, dim_perp = dim
if dim_uv is None:
dim_uv = (dim_rot, max(dim_perp, dim_proj)) # x-y-plane
dim_v, dim_u = dim_uv # y, x
assert dim_v >= dim_rot and dim_u >= dim_perp, 'Projected dimensions are too small!'
# Creating coordinate list of all voxels (for one slice):
voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-x-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj / 2., dim_perp / 2.)
positions = self._get_position(voxels, center, tilt, dim_u)
# Calculate weight-matrix:
r = 1 / np.sqrt(np.pi) # radius of the voxel circle
rho = 0.5 / r
row = []
col = []
data = []
# One slice:
disable = not verbose
for i, voxel in enumerate(tqdm(voxels, disable=disable, leave=False,
desc='Set up projector')):
impacts = self._get_impact(positions[i], r, dim_u) # impact along projected x-axis
voxel_index = voxel[0] * dim_perp * dim_rot + voxel[1] # 0: z, 1: x
for impact in impacts:
impact_index = impact + (dim_v - dim_rot) // 2 * dim_u
distance = np.abs(impact + 0.5 - positions[i])
delta = distance / r
col.append(voxel_index)
row.append(impact_index)
data.append(self._get_weight(delta, rho))
# All other slices (along y):
data = np.tile(data, dim_rot)
columns = np.tile(col, dim_rot)
rows = np.tile(row, dim_rot)
addition = np.repeat(np.arange(dim_rot), len(row))
columns += addition * dim_perp
rows += addition * dim_u
# Calculate weight matrix and coefficients for jacobi matrix:
shape = (np.prod(dim_uv), np.prod(dim))
weight = csr_matrix(coo_matrix((data, (rows, columns)), shape=shape))
coeff = [[np.cos(tilt), 0, np.sin(tilt)], [0, 1, 0]]
super().__init__(dim, dim_uv, weight, coeff)
self._log.debug('Created ' + str(self))
@staticmethod
def _get_position(points, center, tilt, size):
point_vecs = np.asarray(points) + 0.5 - np.asarray(center) # vectors pointing to points
direc_vec = np.array((np.cos(tilt), -np.sin(tilt))) # vector pointing along projection
distances = np.cross(direc_vec, point_vecs) # here (special case): divisor is one!
distances += size / 2. # Shift to the center of the projection
return distances
@staticmethod
def _get_impact(pos, r, size):
return [x for x in np.arange(np.floor(pos - r), np.floor(pos + r) + 1, dtype=int)
if 0 <= x < size]
@staticmethod
def _get_weight(delta, rho): # use circles to represent the voxels
lo, up = delta - rho, delta + rho
# Upper boundary:
if up >= 1:
w_up = 0.5
else:
w_up = (up * np.sqrt(1 - up ** 2) + np.arctan(up / np.sqrt(1 - up ** 2))) / pi
# Lower boundary:
if lo <= -1:
w_lo = -0.5
else:
w_lo = (lo * np.sqrt(1 - lo ** 2) + np.arctan(lo / np.sqrt(1 - lo ** 2))) / pi
return w_up - w_lo
def get_info(self, verbose=False):
"""Get specific information about the projector as a string.
Parameters
----------
verbose: boolean, optional
If this is true, the text looks prettier (maybe using latex). Default is False for the
use in file names and such.
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
"""
if verbose:
return u'y-tilt: $\phi = {:d}$°'.format(int(np.round(self.tilt * 180 / pi)))
else:
return u'ytilt_phi={:d}°'.format(int(np.round(self.tilt * 180 / pi)))
class SimpleProjector(Projector):
"""Class representing a projection function along one of the major axes.
The :class:`~.SimpleProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
axis : {'z', 'y', 'x'}, optional
Main axis along which the magnetic distribution is projected (given as a string). Defaults
to the z-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set it uses the 3D default dimensions.
"""
_log = logging.getLogger(__name__ + '.SimpleProjector')
AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (2, 1, 0)} # (0:z, 1:y, 2:x) -> (proj, v, u)
# coordinate switch for 'x': u, v --> z, y (not y, z!)!
def __init__(self, dim, axis='z', dim_uv=None, verbose=False):
self._log.debug('Calling __init__')
assert axis in {'z', 'y', 'x'}, 'Projection axis has to be x, y or z (given as a string)!'
self.axis = axis
proj, v, u = self.AXIS_DICT[axis]
dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
dim_z, dim_y, dim_x = dim
size_2d = dim_u * dim_v
size_3d = np.prod(dim)
data = np.repeat(1, size_3d) # size_3d ones in the matrix (each voxel is projected)
indptr = np.arange(0, size_3d + 1, dim_proj) # each row has dim_proj 1-entries
if axis == 'z':
self._log.debug('Projecting along the z-axis')
coeff = [[1, 0, 0], [0, 1, 0]]
indices = np.array([np.arange(row, size_3d, size_2d)
for row in range(size_2d)]).reshape(-1)
elif axis == 'y':
self._log.debug('Projection along the y-axis')
coeff = [[1, 0, 0], [0, 0, 1]]
indices = np.array(
[np.arange(row % dim_x, dim_x * dim_y, dim_x) + row // dim_x * dim_x * dim_y
for row in range(size_2d)]).reshape(-1)
elif axis == 'x':
self._log.debug('Projection along the x-axis')
# TODO: is coordinate switch really necessary? Better other way???
coeff = [[0, 0, 1], [0, 1, 0]] # Caution, coordinate switch: u, v --> z, y (not y, z!)
indices = np.array(
[np.arange(dim_x) + (row % dim_z) * dim_x * dim_y + row // dim_z * dim_x
for row in range(size_2d)]).reshape(-1)
else:
raise ValueError('{} is not a valid axis parameter (use x, y or z)!'.format(axis))
if dim_uv is not None:
indptr = list(indptr) # convert to use insert() and append()
# Calculate padding:
d_v = (np.floor((dim_uv[0] - dim_v) / 2).astype(int),
np.ceil((dim_uv[0] - dim_v) / 2).astype(int))
d_u = (np.floor((dim_uv[1] - dim_u) / 2).astype(int),
np.ceil((dim_uv[1] - dim_u) / 2).astype(int))
indptr.extend([indptr[-1]] * d_v[1] * dim_uv[1]) # add empty lines at the end
for i in np.arange(dim_v, 0, -1): # all slices in between
up, lo = i * dim_u, (i - 1) * dim_u # upper / lower slice end
indptr[up:up] = [indptr[up]] * d_u[1] # end of the slice
indptr[lo:lo] = [indptr[lo]] * d_u[0] # start of the slice
indptr = [0] * d_v[0] * dim_uv[1] + indptr # insert empty rows at the beginning
else: # Make sure dim_uv is defined (used for the assertion)
dim_uv = dim_v, dim_u
assert dim_uv[0] >= dim_v and dim_uv[1] >= dim_u, 'Projected dimensions are too small!'
# Create weight-matrix:
shape = (np.prod(dim_uv), np.prod(dim))
weight = csr_matrix((data, indices, indptr), shape=shape)
super().__init__(dim, dim_uv, weight, coeff)
self._log.debug('Created ' + str(self))
def get_info(self, verbose=False):
"""Get specific information about the projector as a string.
Parameters
----------
verbose: boolean, optional
If this is true, the text looks prettier (maybe using latex). Default is False for the
use in file names and such.
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
"""
if verbose:
return 'projected along {}-axis'.format(self.axis)
else:
return '{}axis'.format(self.axis)
pyramid/ramp.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.Ramp` class which implements polynomial phase ramps."""
import numpy as np
from pyramid.phasemap import PhaseMap
__all__ = ['Ramp']
class Ramp(object):
"""Class representing a polynomial phase ramp.
Sometimes additional phase ramps occur in phase maps which do not stem from a magnetization
distribution inside the FOV. This class allows the construction (and via the derivative
functions also the reconstruction) of a polynomial ramp. This class is generally constructed
within the ForwardModel and can be retrieved as its attribute if ramp information should be
accessed.
Attributes
----------
data_set : :class:`~dataset.DataSet`
:class:`~dataset.DataSet` object, which stores all required information calculation.
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).
deg_of_freedom : int
Number of degrees of freedom. This is calculated to ``1 + 2 * order``. There is just one
degree of freedom for a ramp of order zero (offset), every higher order contributes two
degrees of freedom.
param_cache : :class:`numpy.ndarray` (N=2)
Parameter cache which is used to store the polynomial coefficients. Higher coefficients
(one for each degree of freedom) are saved along the first axis, values for different
images along the second axis.
n : int
Size of the input space. Coincides with the numer of entries in `param_cache` and
calculates to ``deg_of_freedom * data_set.count``.
Notes
-----
After a reconstruction the relevant polynomial ramp information is stored in the
`param_cache`. If a phasemap with index `i` in the DataSet should be corrected use:
.. code-block:: python
phasemap -= ramp(i=0, dof_list=[0, 1, 2])
The optional parameter `dof_list` can be used to specify a list of degrees of freedom which
should be used for the ramp (e.g. `[0]` will just apply the offset, `[0, 1, 2]` will apply
the offset and linear ramps in both directions).
Fitting polynoms of higher orders than `order = 1` is possible but not recommended, because
features which stem from the magnetization could be covered by the polynom, decreasing the
phase contribution of the magnetization distribution, leading to a false retrieval.
"""
def __init__(self, data_set, order=None):
assert order is None or (isinstance(order, int) and order >= 0), \
'Order has to be None or a positive integer!'
self.order = order
self.a = data_set.a
self.count = data_set.count
self.dimensions = [projector.dim_uv for projector in data_set.projectors]
self.hook_points = data_set.hook_points
self.deg_of_freedom = (1 + 2 * self.order) if self.order is not None else 0
self.param_cache = np.zeros((self.deg_of_freedom, self.count))
self.n = self.deg_of_freedom * self.count # 0 if order is None
def __call__(self, index, dof_list=None):
if self.order is None: # Do nothing if order is None!
return 0
else:
if dof_list is None: # if no specific list is supplied!
dof_list = range(self.deg_of_freedom) # use all available degrees of freedom
dim_uv = self.dimensions[index]
phase_ramp = np.zeros(dim_uv)
# Iterate over all degrees of freedom:
for dof in dof_list:
# Add the contribution of the current degree of freedom:
phase_ramp += (self.param_cache[dof][index] *
self.create_poly_mesh(self.a, dof, dim_uv))
return PhaseMap(self.a, phase_ramp, mask=np.zeros(dim_uv, dtype=np.bool))
def jac_dot(self, index):
"""Calculate the product of the Jacobi matrix .
Parameters
----------
index : int
Index of the phasemap from the `dataset` for which the phase ramp is calculated.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the input
`vector`. Just the ramp contribution is calculated!
"""
if self.order is None: # Do nothing if order is None!
return 0
else:
dim_uv = self.dimensions[index]
phase_ramp = np.zeros(dim_uv)
# Iterate over all degrees of freedom:
for dof in range(self.deg_of_freedom):
# Add the contribution of the current degree of freedom:
phase_ramp += (self.param_cache[dof][index] *
self.create_poly_mesh(self.a, dof, dim_uv))
return np.ravel(phase_ramp)
def jac_T_dot(self, vector):
"""'Calculate the transposed ramp parameters from a given `vector`.
Parameters
----------
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)
Transposed ramp parameters.
"""
result = []
hp = self.hook_points
# Iterate over all degrees of freedom:
for dof in range(self.deg_of_freedom):
# Iterate over all projectors:
for i, dim_uv in enumerate(self.dimensions):
sub_vec = vector[hp[i]:hp[i + 1]]
poly_mesh = self.create_poly_mesh(self.a, dof, dim_uv)
# Transposed ramp parameters: summed product of the vector with the poly-mesh:
result.append(np.sum(sub_vec * np.ravel(poly_mesh)))
return result
def extract_ramp_params(self, x):
"""Extract the ramp parameters of an input vector and return the rest.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Input vector which consists of the vectorised magnetization distribution and the ramp
parameters at the end which will be extracted.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Inpput vector without the extracted ramp parameters.
Notes
-----
This method should always be used before a vector `x` is processed if it is known that
ramp parameters are present so that other functions do not have to bother with them
and the :class:`.~ramp` already knows all important parameters for its own functions.
"""
if self.order is not None: # Do nothing if order is None!
# Split off ramp parameters and fill cache:
x, ramp_params = np.split(x, [-self.n])
self.param_cache = ramp_params.reshape((self.deg_of_freedom, self.count))
return x
@classmethod
def create_poly_mesh(cls, a, deg_of_freedom, dim_uv):
"""Create a polynomial mesh for the ramp calculation for a specific degree of freedom.
Parameters
----------
a : float
Grid spacing which should be used for the ramp.
deg_of_freedom : int
Current degree of freedom for which the mesh should be created. 0 corresponds to a
simple offset, 1 corresponds to a linear ramp in u-direction, 2 to a linear ramp in
v-direction and so on.
dim_uv : tuple (N=2)
Dimensions of the 2D mesh that should be created.
Returns
-------
result_mesh : :class:`~numpy.ndarray` (N=2)
Polynomial mesh that was created and can be used for further calculations.
"""
# Determine if u-direction (u_or_v == 1) or v-direction (u_or_v == 0)!
u_or_v = (deg_of_freedom - 1) % 2
# Determine polynomial order:
order = (deg_of_freedom + 1) // 2
# Return polynomial mesh:
return (np.indices(dim_uv)[u_or_v] * a) ** order
@classmethod
def create_ramp(cls, a, dim_uv, params):
"""Class method to create an arbitrary polynomial ramp.
Parameters
----------
a : float
Grid spacing which should be used for the ramp.
dim_uv : tuple (N=2)
Dimensions of the 2D mesh that should be created.
params : list
List of ramp parameters. The first entry corresponds to a simple offset, the second
and third correspond to a linear ramp in v- and u-direction, respectively and so on.
Returns
-------
phase_ramp : :class:`~pyramid.phasemap.PhaseMap`
The phase ramp as a :class:`~pyramid.phasemap.PhaseMap` object.
"""
phase_ramp = np.zeros(dim_uv)
dof_list = range(len(params))
for dof in dof_list:
phase_ramp += params[dof] * cls.create_poly_mesh(a, dof, dim_uv)
# Return the phase ramp as a PhaseMap with empty (!) mask:
return PhaseMap(a, phase_ramp, mask=np.zeros(dim_uv, dtype=np.bool))
pyramid/reconstruction.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Reconstruct magnetic distributions from given phasemaps.
This module reconstructs 3-dimensional magnetic distributions (as
:class:`~pyramid.magdata.VectorData` objects) from a given set of phase maps (represented by
:class:`~pyramid.phasemap.PhaseMap` objects) by using several model based reconstruction algorithms
which use the forward model provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper`
and a priori knowledge of the distribution.
"""
import logging
import numpy as np
from pyramid.fielddata import VectorData
__all__ = ['optimize_linear', 'optimize_nonlin', 'optimize_splitbregman']
_log = logging.getLogger(__name__)
def optimize_linear(costfunction, mag_0=None, ramp_0=None, max_iter=None, verbose=False):
"""Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
Blazingly fast for l2-based cost functions.
Parameters
----------
costfunction : :class:`~.Costfunction`
A :class:`~.Costfunction` object which implements a specified forward model and
regularisator which is minimized in the optimization process.
mag_0: :class:`~.VectorData`
The starting magnetisation distribution used for the reconstruction. A zero vector will be
used if no VectorData object is specified.
mag_0: :class:`~.Ramp`
The starting ramp for the reconstruction. A zero vector will be
used if no Ramp object is specified.
max_iter : int, optional
The maximum number of iterations for the opimization.
verbose: bool, optional
If set to True, information like a progressbar is displayed during reconstruction.
The default is False.
Returns
-------
magdata : :class:`~pyramid.fielddata.VectorData`
The reconstructed magnetic distribution as a :class:`~.VectorData` object.
"""
import jutil.cg as jcg
from jutil.taketime import TakeTime
_log.debug('Calling optimize_linear')
_log.info('Cost before optimization: {:.3e}'.format(costfunction(np.zeros(costfunction.n))))
data_set = costfunction.fwd_model.data_set
# Get starting distribution vector x_0:
x_0 = np.empty(costfunction.n)
if mag_0 is not None:
costfunction.fwd_model.magdata = mag_0
x_0[:data_set.n] = costfunction.fwd_model.magdata.get_vector(mask=data_set.mask)
if ramp_0 is not None:
ramp_vec = ramp_0.param_cache.ravel()
else:
ramp_vec = np.zeros_like(costfunction.fwd_model.ramp.n)
x_0[data_set.n:] = ramp_vec
# Minimize:
with TakeTime('reconstruction time'):
x_opt = jcg.conj_grad_minimize(costfunction, x_0=x_0, max_iter=max_iter, verbose=verbose).x
_log.info('Cost after optimization: {:.3e}'.format(costfunction(x_opt)))
# Cut ramp parameters if necessary (this also saves the final parameters in the ramp class!):
x_opt = costfunction.fwd_model.ramp.extract_ramp_params(x_opt)
# Create and return fitting VectorData object:
mag_opt = VectorData(data_set.a, np.zeros((3,) + data_set.dim))
mag_opt.set_vector(x_opt, data_set.mask)
return mag_opt
def optimize_nonlin(costfunction, first_guess=None):
"""Reconstruct a three-dimensional magnetic distribution from given phase maps via
steepest descent method. This is slow, but works best for non l2-regularisators.
Parameters
----------
costfunction : :class:`~.Costfunction`
A :class:`~.Costfunction` object which implements a specified forward model and
regularisator which is minimized in the optimization process.
first_guess : :class:`~pyramid.fielddata.VectorData`
magnetization to start the non-linear iteration with.
Returns
-------
magdata : :class:`~pyramid.fielddata.VectorData`
The reconstructed magnetic distribution as a :class:`~.VectorData` object.
"""
import jutil.minimizer as jmin
import jutil.norms as jnorms
_log.debug('Calling optimize_nonlin')
data_set = costfunction.fwd_model.data_set
if first_guess is None:
first_guess = VectorData(data_set.a, np.zeros((3,) + data_set.dim))
x_0 = first_guess.get_vector(data_set.mask)
assert len(x_0) == costfunction.n, (len(x_0), costfunction.m, costfunction.n)
p = costfunction.regularisator.p
q = 1. / (1. - (1. / p))
lq = jnorms.LPPow(q, 1e-20)
def _preconditioner(_, direc):
direc_p = direc / abs(direc).max()
direc_p = 10 * (1. / q) * lq.jac(direc_p)
return direc_p
# This Method is semi-best for Lp type problems. Takes forever, though
_log.info('Cost before optimization: {}'.format(costfunction(np.zeros(costfunction.n))))
result = jmin.minimize(
costfunction, x_0,
method="SteepestDescent",
options={"preconditioner": _preconditioner},
tol={"max_iteration": 10000})
x_opt = result.x
_log.info('Cost after optimization: {}'.format(costfunction(x_opt)))
mag_opt = VectorData(data_set.a, np.zeros((3,) + data_set.dim))
mag_opt.set_vector(x_opt, data_set.mask)
return mag_opt
def optimize_splitbregman(costfunction, weight, lam, mu):
"""
Reconstructs magnet distribution from phase image measurements using a split bregman
algorithm with a dedicated TV-l1 norm. Very dedicated, frickle, brittle, and difficult
to get to work, but fastest option available if it works.
Seems to work for some 2D examples with weight=lam=1 and mu in [1, .., 1e4].
Parameters
----------
costfunction : :class:`~.Costfunction`
A :class:`~.Costfunction` object which implements a specified forward model and
regularisator which is minimized in the optimization process.
weight : float
Obscure split bregman parameter
lam : float
Cryptic split bregman parameter
mu : float
flabberghasting split bregman paramter
Returns
-------
magdata : :class:`~pyramid.fielddata.VectorData`
The reconstructed magnetic distribution as a :class:`~.VectorData` object.
"""
import jutil.splitbregman as jsb
import jutil.operator as joperator
import jutil.diff as jdiff
_log.debug('Calling optimize_splitbregman')
# regularisator is actually not necessary, but this makes the cost
# function to that which is supposedly optimized by split bregman.
# Thus cost can be used to verify convergence
fwd_model = costfunction.fwd_model
data_set = fwd_model.data_set
A = joperator.Function(
(costfunction.m, costfunction.n),
lambda x: fwd_model.jac_dot(None, x),
FT=lambda x: fwd_model.jac_T_dot(None, x))
D = joperator.VStack([
jdiff.get_diff_operator(data_set.mask, 0, 3),
jdiff.get_diff_operator(data_set.mask, 1, 3)])
y = np.asarray(costfunction.y, dtype=np.double)
x_opt = jsb.split_bregman_2d(
A, D, y,
weight=weight, mu=mu, lambd=lam, max_iter=1000)
mag_opt = VectorData(data_set.a, np.zeros((3,) + data_set.dim))
mag_opt.set_vector(x_opt, data_set.mask)
return mag_opt
pyramid/regularisator.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides the :class:`~.Regularisator` class which represents a regularisation term
which adds additional constraints to a costfunction to minimize."""
import abc
import logging
import numpy as np
from scipy import sparse
import jutil.diff as jdiff
import jutil.norms as jnorm
__all__ = ['NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator',
'ComboRegularisator']
class Regularisator(object, metaclass=abc.ABCMeta):
"""Class for providing a regularisation term which implements additional constraints.
Represents a certain constraint for the 3D magnetization distribution whose cost is to minimize
in addition to the derivation from the 2D phase maps. Important is the used `norm` and the
regularisation parameter `lam` (lambda) which determines the weighting between the two cost
parts (measurements and regularisation). Additional parameters at the end of the input
vector, which are not relevant for the regularisation can be discarded by specifying the
number in `add_params`.
Attributes
----------
norm : :class:`~jutil.norm.WeightedNorm`
Norm, which is used to determine the cost of the regularisation term.
lam : float
Regularisation parameter determining the weighting between measurements and regularisation.
add_params : int
Number of additional parameters which are not used in the regularisation. Used to cut
the input vector into the appropriate size.
"""
_log = logging.getLogger(__name__ + '.Regularisator')
@abc.abstractmethod
def __init__(self, norm, lam, add_params=0):
self._log.debug('Calling __init__')
self.norm = norm
self.lam = lam
self.add_params = add_params
if self.add_params > 0:
self.slice = slice(-add_params)
else:
self.slice = slice(None)
self._log.debug('Created ' + str(self))
def __call__(self, x):
self._log.debug('Calling __call__')
return self.lam * self.norm(x[self.slice])
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(norm=%r, lam=%r, add_params=%r)' % (self.__class__, self.norm, self.lam,
self.add_params)
def __str__(self):
self._log.debug('Calling __str__')
return 'Regularisator(norm=%s, lam=%s, add_params=%s)' % (self.norm, self.lam,
self.add_params)
def jac(self, x):
"""Calculate the derivative of the regularisation term for a given magnetic 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.
"""
result = np.zeros_like(x)
result[self.slice] = self.lam * self.norm.jac(x[self.slice])
return result
def hess_dot(self, x, vector):
"""Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
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.
"""
result = np.zeros_like(vector)
result[self.slice] = self.lam * self.norm.hess_dot(x, vector[self.slice])
return result
def hess_diag(self, x):
""" Return the diagonal of the Hessian.
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 in the
computation). It is implemented for the case that in the future nonlinear problems
have to be solved.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Diagonal of the Hessian matrix.
"""
self._log.debug('Calling hess_diag')
result = np.zeros_like(x)
result[self.slice] = self.lam * self.norm.hess_diag(x[self.slice])
return result
class ComboRegularisator(Regularisator):
"""Class for providing a regularisation term which combines several regularisators.
If more than one regularisation should be utilized, this class can be use. It is given a list
of :class:`~.Regularisator` objects. The input will be forwarded to each of them and the
results are summed up and returned.
Attributes
----------
reg_list: :class:`~.Regularisator`
A list of regularisator objects to whom the input is passed on.
"""
def __init__(self, reg_list):
self._log.debug('Calling __init__')
self.reg_list = reg_list
super().__init__(norm=None, lam=None)
self._log.debug('Created ' + str(self))
def __call__(self, x):
self._log.debug('Calling __call__')
return np.sum([self.reg_list[i](x) for i in range(len(self.reg_list))], axis=0)
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(reg_list=%r)' % (self.__class__, self.reg_list)
def __str__(self):
self._log.debug('Calling __str__')
return 'ComboRegularisator(reg_list=%s)' % self.reg_list
def jac(self, x):
"""Calculate the derivative of the regularisation term for a given magnetic 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.
"""
return np.sum([self.reg_list[i].jac(x) for i in range(len(self.reg_list))], axis=0)
def hess_dot(self, x, vector):
"""Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
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.
"""
return np.sum([self.reg_list[i].hess_dot(x, vector) for i in range(len(self.reg_list))],
axis=0)
def hess_diag(self, x):
""" Return the diagonal of the Hessian.
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 in the
computation). It is implemented for the case that in the future nonlinear problems
have to be solved.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Diagonal of the Hessian matrix.
"""
self._log.debug('Calling hess_diag')
return np.sum([self.reg_list[i].hess_diag(x) for i in range(len(self.reg_list))], axis=0)
class NoneRegularisator(Regularisator):
"""Placeholder class if no regularization is used.
This class is instantiated in the :class:`~pyramid.costfunction.Costfunction`, which means
no regularisation is used. All associated functions return appropriate zero-values.
Attributes
----------
norm: None
No regularization is used, thus also no norm.
lam: 0
Not used.
"""
_log = logging.getLogger(__name__ + '.NoneRegularisator')
def __init__(self):
self._log.debug('Calling __init__')
self.norm = None
self.lam = 0
self.add_params = None
super().__init__(norm=None, lam=None)
self._log.debug('Created ' + str(self))
def __call__(self, x):
self._log.debug('Calling __call__')
return 0
def jac(self, x):
"""Calculate the derivative of the regularisation term for a given magnetic 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.
"""
return np.zeros_like(x)
def hess_dot(self, x, vector):
"""Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
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 in the
computation). It is implemented for the case that in the future nonlinear problems
have to be solved.
vector : a :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 np.zeros_like(vector)
def hess_diag(self, x):
""" Return the diagonal of the Hessian.
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.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Diagonal of the Hessian matrix.
"""
self._log.debug('Calling hess_diag')
return np.zeros_like(x)
class ZeroOrderRegularisator(Regularisator):
"""Class for providing a regularisation term which implements Lp norm minimization.
The constraint this class represents is the minimization of the Lp norm for the 3D
magnetization distribution. Important is the regularisation parameter `lam` (lambda) which
determines the weighting between the two cost parts (measurements and regularisation).
Attributes
----------
lam: float
Regularisation parameter determining the weighting between measurements and regularisation.
p: int, optional
Order of the norm (default: 2, which means a standard L2-norm).
add_params : int
Number of additional parameters which are not used in the regularisation. Used to cut
the input vector into the appropriate size.
"""
_log = logging.getLogger(__name__ + '.ZeroOrderRegularisator')
def __init__(self, _=None, lam=1E-4, p=2, add_params=0):
self._log.debug('Calling __init__')
self.p = p
if p == 2:
norm = jnorm.L2Square()
else:
norm = jnorm.LPPow(p, 1e-12)
super().__init__(norm, lam, add_params)
self._log.debug('Created ' + str(self))
class FirstOrderRegularisator(Regularisator):
"""Class for providing a regularisation term which implements derivation minimization.
The constraint this class represents is the minimization of the first order derivative of the
3D magnetization distribution using a Lp norm. Important is the regularisation parameter `lam`
(lambda) which determines the weighting between the two cost parts (measurements and
regularisation).
Attributes
----------
mask: :class:`~numpy.ndarray` (N=3)
A boolean mask which defines the magnetized volume in 3D.
lam: float
Regularisation parameter determining the weighting between measurements and regularisation.
p: int, optional
Order of the norm (default: 2, which means a standard L2-norm).
add_params : int
Number of additional parameters which are not used in the regularisation. Used to cut
the input vector into the appropriate size.
"""
def __init__(self, mask, lam=1E-4, p=2, add_params=0):
self.p = p
D0 = jdiff.get_diff_operator(mask, 0, 3)
D1 = jdiff.get_diff_operator(mask, 1, 3)
D2 = jdiff.get_diff_operator(mask, 2, 3)
D = sparse.vstack([D0, D1, D2])
if p == 2:
norm = jnorm.WeightedL2Square(D)
else:
norm = jnorm.WeightedTV(jnorm.LPPow(p, 1e-12), D, [D0.shape[0], D.shape[0]])
super().__init__(norm, lam, add_params)
self._log.debug('Created ' + str(self))
pyramid/tests/.DS_Store deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_analytic.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
"""Testcase for the analytic module."""
import os
import unittest
import numpy as np
from numpy import pi
from numpy.testing import assert_allclose
import pyramid.analytic as an
class TestCaseAnalytic(unittest.TestCase):
"""TestCase for the analytic module."""
path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_analytic/')
dim = (4, 4, 4)
a = 10.0
phi = pi / 4
center = (dim[0] / 2, dim[1] / 2, dim[2] / 2)
radius = dim[2] / 4
def test_phase_mag_slab(self):
"""Test of the phase_mag_slab method."""
width = (self.dim[0] / 2, self.dim[1] / 2, self.dim[2] / 2)
phase = an.phase_mag_slab(self.dim, self.a, self.phi, self.center, width).phase
reference = np.load(os.path.join(self.path, 'ref_phase_slab.npy'))
assert_allclose(phase, reference, atol=1E-10,
err_msg='Unexpected behavior in phase_mag_slab()')
def test_phase_mag_disc(self):
"""Test of the phase_mag_disc method."""
radius = self.dim[2] / 4
height = self.dim[2] / 2
phase = an.phase_mag_disc(self.dim, self.a, self.phi, self.center, radius, height).phase
reference = np.load(os.path.join(self.path, 'ref_phase_disc.npy'))
assert_allclose(phase, reference, atol=1E-10,
err_msg='Unexpected behavior in phase_mag_disc()')
def test_phase_mag_sphere(self):
"""Test of the phase_mag_sphere method."""
radius = self.dim[2] / 4
phase = an.phase_mag_sphere(self.dim, self.a, self.phi, self.center, radius).phase
reference = np.load(os.path.join(self.path, 'ref_phase_sphere.npy'))
assert_allclose(phase, reference, atol=1E-10,
err_msg='Unexpected behavior in phase_mag_sphere()')
def test_phase_mag_vortex(self):
"""Test of the phase_mag_vortex method."""
radius = self.dim[2] / 4
height = self.dim[2] / 2
phase = an.phase_mag_vortex(self.dim, self.a, self.center, radius, height).phase
reference = np.load(os.path.join(self.path, 'ref_phase_vort.npy'))
assert_allclose(phase, reference, atol=1E-10,
err_msg='Unexpected behavior in phase_mag_vortex()')
pyramid/tests/test_analytic/ref_phase_disc.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_analytic/ref_phase_slab.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_analytic/ref_phase_sphere.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_analytic/ref_phase_vort.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_costfunction.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
"""Testcase for the costfunction module"""
import os
import unittest
import numpy as np
from numpy.testing import assert_allclose
from pyramid.costfunction import Costfunction
from pyramid.dataset import DataSet
from pyramid.forwardmodel import ForwardModel
from pyramid.projector import SimpleProjector
from pyramid.regularisator import FirstOrderRegularisator
from pyramid import load_phasemap
class TestCaseCostfunction(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_costfunction')
self.a = 10.
self.dim = (4, 5, 6)
self.mask = np.zeros(self.dim, dtype=bool)
self.mask[1:-1, 1:-1, 1:-1] = True
self.data = DataSet(self.a, self.dim, mask=self.mask)
self.projector = SimpleProjector(self.dim)
self.phasemap = load_phasemap(os.path.join(self.path, 'phasemap_ref.hdf5'))
self.data.append(self.phasemap, self.projector)
self.data.append(self.phasemap, self.projector)
self.reg = FirstOrderRegularisator(self.mask, lam=1E-4)
self.cost = Costfunction(ForwardModel(self.data), self.reg)
def tearDown(self):
self.path = None
self.a = None
self.dim = None
self.mask = None
self.data = None
self.projector = None
self.phasemap = None
self.reg = None
self.cost = None
def test_call(self):
assert_allclose(self.cost(np.ones(self.cost.n)), 0., atol=1E-7,
err_msg='Unexpected behaviour in __call__()!')
zero_vec_cost = np.load(os.path.join(self.path, 'zero_vec_cost.npy'))
assert_allclose(self.cost(np.zeros(self.cost.n)), zero_vec_cost,
err_msg='Unexpected behaviour in __call__()!')
def test_jac(self):
assert_allclose(self.cost.jac(np.ones(self.cost.n)), np.zeros(self.cost.n), atol=1E-7,
err_msg='Unexpected behaviour in jac()!')
jac_vec_ref = np.load(os.path.join(self.path, 'jac_vec_ref.npy'))
assert_allclose(self.cost.jac(np.zeros(self.cost.n)), jac_vec_ref, atol=1E-7,
err_msg='Unexpected behaviour in jac()!')
jac = np.array([self.cost.jac(np.eye(self.cost.n)[:, i]) for i in range(self.cost.n)]).T
jac_ref = np.load(os.path.join(self.path, 'jac_ref.npy'))
assert_allclose(jac, jac_ref, atol=1E-7,
err_msg='Unexpected behaviour in jac()!')
def test_hess_dot(self):
assert_allclose(self.cost.hess_dot(None, np.zeros(self.cost.n)), np.zeros(self.cost.n),
atol=1E-7, err_msg='Unexpected behaviour in jac()!')
hess_vec_ref = -np.load(os.path.join(self.path, 'jac_vec_ref.npy'))
assert_allclose(self.cost.hess_dot(None, np.ones(self.cost.n)), hess_vec_ref, atol=1E-7,
err_msg='Unexpected behaviour in jac()!')
hess = np.array([self.cost.hess_dot(None, np.eye(self.cost.n)[:, i])
for i in range(self.cost.n)]).T
hess_ref = np.load(os.path.join(self.path, 'hess_ref.npy'))
assert_allclose(hess, hess_ref, atol=1E-7,
err_msg='Unexpected behaviour in hess_dot()!')
def test_hess_diag(self):
assert_allclose(self.cost.hess_diag(None), np.ones(self.cost.n),
err_msg='Unexpected behaviour in hess_diag()!')
pyramid/tests/test_costfunction/hess_ref.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_costfunction/jac_ref.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_costfunction/jac_vec_ref.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_costfunction/phasemap_ref.hdf5 deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_costfunction/zero_vec_cost.npy deleted 100644 → 0
View file @ 45d1150a
File deleted
pyramid/tests/test_dataset.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
"""Testcase for the dataset module"""
import os
import unittest
import numpy as np
from numpy.testing import assert_allclose
from pyramid.dataset import DataSet
from pyramid.fielddata import VectorData
from pyramid.phasemap import PhaseMap
from pyramid.projector import SimpleProjector
class TestCaseDataSet(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_dataset')
self.a = 10.
self.dim = (4, 5, 6)
self.mask = np.zeros(self.dim, dtype=bool)
self.mask[1:-1, 1:-1, 1:-1] = True
self.data = DataSet(self.a, self.dim, mask=self.mask)
self.projector = SimpleProjector(self.dim)
self.phasemap = PhaseMap(self.a, np.ones(self.dim[1:3]))
def tearDown(self):
self.path = None
self.a = None
self.dim = None
self.mask = None
self.data = None
self.projector = None
self.phasemap = None
def test_append(self):
self.data.append(self.phasemap, self.projector)
assert self.data.phasemaps[0] == self.phasemap, 'Phase map not correctly assigned!'
assert self.data.projectors[0] == self.projector, 'Projector not correctly assigned!'
def test_create_phasemaps(self):
self.data.append(PhaseMap(self.a, np.zeros(self.projector.dim_uv)), self.projector)
magdata = VectorData(self.a, np.ones((3,) + self.dim))
phasemaps = self.data.create_phasemaps(magdata)
phase_vec = phasemaps[0].phase_vec
phase_vec_ref = np.load(os.path.join(self.path, 'phase_vec_ref.npy'))
assert_allclose(phase_vec, phase_vec_ref, atol=1E-6,
err_msg='Unexpected behaviour in create_phasemaps()!')
def test_set_Se_inv_block_diag(self):
self.data.append(self.phasemap, self.projector)
self.data.append(self.phasemap, self.projector)
cov = np.diag(np.ones(np.prod(self.phasemap.dim_uv)))
self.data.set_Se_inv_block_diag([cov, cov])
assert self.data.Se_inv.shape == (self.data.m, self.data.m), \
'Unexpected behaviour in set_Se_inv_block_diag()!'
assert self.data.Se_inv.diagonal().sum() == self.data.m, \
'Unexpected behaviour in set_Se_inv_block_diag()!'
def test_set_Se_inv_diag_with_conf(self):
self.data.append(self.phasemap, self.projector)
self.data.append(self.phasemap, self.projector)
confidence = self.mask[1, ...]
self.data.set_Se_inv_diag_with_conf([confidence, confidence])
assert self.data.Se_inv.shape == (self.data.m, self.data.m), \
'Unexpected behaviour in set_Se_inv_diag_with_masks()!'
assert self.data.Se_inv.diagonal().sum() == 2 * confidence.sum(), \
'Unexpected behaviour in set_Se_inv_diag_with_masks()!'
def test_set_3d_mask(self):
projector_z = SimpleProjector(self.dim, axis='z')
projector_y = SimpleProjector(self.dim, axis='y')
projector_x = SimpleProjector(self.dim, axis='x')
mask_z = np.zeros(projector_z.dim_uv, dtype=bool)
mask_y = np.zeros(projector_y.dim_uv, dtype=bool)
mask_x = np.zeros(projector_x.dim_uv, dtype=bool)
mask_z[1:-1, 1:-1] = True
mask_y[1:-1, 1:-1] = True
mask_x[1:-1, 1:-1] = True
phasemap_z = PhaseMap(self.a, np.zeros(projector_z.dim_uv), mask_z)
phasemap_y = PhaseMap(self.a, np.zeros(projector_y.dim_uv), mask_y)
phasemap_x = PhaseMap(self.a, np.zeros(projector_x.dim_uv), mask_x)
self.data.append(phasemap_z, projector_z)
self.data.append(phasemap_y, projector_y)
self.data.append(phasemap_x, projector_x)
self.data.set_3d_mask()
mask_ref = np.zeros(self.dim, dtype=bool)
mask_ref[1:-1, 1:-1, 1:-1] = True
np.testing.assert_equal(self.data.mask, mask_ref,
err_msg='Unexpected behaviour in set_3d_mask')