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 56 additions and 4971 deletions
pyramid/phasemap.py deleted 100644 → 0
View file @ f6abe2ad
# -*- 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_lim = phase
# Clip non-trustworthy regions for the limit calculation:
if show_conf:
phase_trust = np.where(self.confidence > 0.9, phase_lim, np.nan)
phase_min, phase_max = np.nanmin(phase_trust), np.nanmax(phase_trust)
phase_lim = np.clip(phase_lim, phase_min, phase_max)
# Cut outlier beyond a certain sigma-margin:
if sigma_clip is not None:
outlier = np.abs(phase_lim - np.mean(phase_lim)) < sigma_clip * np.std(phase_lim)
phase_sigma = np.where(outlier, phase_lim, np.nan)
phase_min, phase_max = np.nanmin(phase_sigma), np.nanmax(phase_sigma)
phase_lim = np.clip(phase_lim, phase_min, phase_max)
# Calculate the limits if necessary (zero has to be present!):
if vmin is None:
vmin = np.min(phase_lim)
if vmax is None:
vmax = np.max(phase_lim)
# 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, sigma_clip=2,
**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:
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 @ f6abe2ad
# -*- 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):
"""Class representing a phase mapping strategy for the electrostatic charge contribution.
The :class:`~.PhaseMapperCharge` class represents a phase mapping strategy for the electrostatic
charge contribution to the electron phase shift which results e.g. from the charges 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
----------
kernelcharge : :class:`~pyramid.KernelCharge`
Convolution kernel, representing the phase contribution of one single charge pixel.
m: int
Size of the image space.
n: int
Size of the input space.
c: :class:`~numpy.ndarray` (N=3)
The charge distribution
c_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the charge distribution.
phase_fft: :class:`~numpy.ndarray` (N=3)
The real FFT of the phase from the given charge distribution c.
"""
_log = logging.getLogger(__name__ + '.PhaseMapperCharge')
def __init__(self, kernelcharge):
self._log.debug('Calling __init__')
self.kernelcharge = kernelcharge
self.m = np.prod(kernelcharge.dim_uv)
self.n = self.m
self.c = np.zeros(kernelcharge.dim_pad, dtype=kernelcharge.kc.dtype)
self.phase_adj = np.zeros(kernelcharge.dim_pad, dtype=kernelcharge.kc.dtype)
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(kernelcharge=%r)' % (self.__class__, self.kernelcharge)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMapperCharge(kernelcharge=%s)' % self.kernelcharge
def __call__(self, elecdata):
assert isinstance(elecdata, ScalarData), 'Only ScalarData objects can be mapped!'
assert elecdata.a == self.kernelcharge.a, 'Grid spacing has to match!'
assert elecdata.dim[0] == 1, 'Charge distribution must be 2-dimensional!'
assert elecdata.dim[1:3] == self.kernelcharge.dim_uv, 'Dimensions do not match!'
# Process input parameters:
self.c[self.kernelcharge.slice_c] = elecdata.field[0, ...]
return PhaseMap(elecdata.a, self._convolve())
def _convolve(self):
# Fourier transform of the projected charge distribution:
self.c_fft = fft.rfftn(self.c)
# Convolve the charge distribution with the kernel in Fourier space:
self.phase_fft = self.c_fft * self.kernelcharge.kc_fft
# Return the result:
return fft.irfftn(self.phase_fft)[self.kernelcharge.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 charge distribution of every pixel (row-wise).
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitly calculated) with the vector.
"""
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
self.c[self.kernelcharge.slice_c] = np.reshape(vector, self.kernelcharge.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 vector
phasemap.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitly calculated) with
the vector, which has ``N**2`` entries like a 2D charge projection.
"""
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
self.phase_adj[self.kernelcharge.slice_phase] = vector.reshape(self.kernelcharge.dim_uv)
phase_adj_fft = fft.irfft2_adj(self.phase_adj)
kc_adj_fft = phase_adj_fft * np.conj(self.kernelcharge.kc_fft)
kc_adj = fft.rfft2_adj(kc_adj_fft)[self.kernelcharge.slice_c]
result = kc_adj.ravel()
return result
pyramid/plottools.py deleted 100644 → 0
View file @ f6abe2ad
# -*- 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
# TODO: Replace by matplotlib styles!
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!
# TODO: set 'font.size' = FONTSIZE! Does this change everything relative to this value???
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 = None#axis.get_facecolor() TODO: Activate for matplotlib 2.0!
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 @ f6abe2ad
# -*- 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 @ f6abe2ad
# -*- 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 @ f6abe2ad
# -*- 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, ScalarData
__all__ = ['optimize_linear', 'optimize_linear_charge', '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 optimization 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.
ramp_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 optimization.
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_linear_charge(costfunction, charge_0=None, ramp_0=None, max_iter=None, verbose=False):
"""Reconstruct a three-dimensional charge distribution from given phase maps via the
conjugate gradient optimization 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.
charge_0: :class:`~.VectorData`
The starting magnetisation distribution used for the reconstruction. A zero vector will be
used if no VectorData object is specified.
ramp_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 optimization.
verbose: bool, optional
If set to True, information like a progressbar is displayed during reconstruction.
The default is False.
Returns
-------
elecdata : :class:`~pyramid.fielddata.ScalarData`
The reconstructed charge distribution as a :class:`~.ScalarData` object.
"""
import jutil.cg as jcg
from jutil.taketime import TakeTime
_log.debug('Calling optimize_linear_charge')
_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 charge_0 is not None:
costfunction.fwd_model.elecdata = charge_0
x_0[:data_set.n] = costfunction.fwd_model.elecdata.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 ScalarData object:
charge_opt = ScalarData(data_set.a, np.zeros(data_set.dim))
charge_opt.set_vector(x_opt, data_set.mask)
return charge_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 @ f6abe2ad
# -*- 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/utils/__init__.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing Pyramid utility functions."""
from .pm import pm
from .reconstruction_2d_from_phasemap import reconstruction_2d_from_phasemap
from .reconstruction_3d_from_magdata import reconstruction_3d_from_magdata
#from .phasemap_creator import gui_phasemap_creator
#from .mag_slicer import gui_mag_slicer
__all__ = ['pm', 'reconstruction_2d_from_phasemap', 'reconstruction_3d_from_magdata']#,
#'gui_phasemap_creator', 'gui_mag_slicer']
pyramid/utils/mag_slicer.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
# Form implementation generated from reading ui file 'mag_slicer.ui'
#
# Created: Sun Aug 31 20:39:52 2014
# by: PyQt4 UI code generator 4.9.6
#
# WARNING! All changes made in this file will be lost!
"""GUI for slicing 3D magnetization distributions."""
import logging
import os
import sys
from matplotlib.figure import Figure
from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
from ..projector import SimpleProjector
from ..kernel import Kernel
from ..phasemapper import PhaseMapperRDFC
from ..file_io.io_vectordata import load_vectordata
try:
from PyQt5 import QtGui, QtCore
from PyQt5.uic import loadUiType
except ImportError:
from PyQt4 import QtGui, QtCore
from PyQt4.uic import loadUiType
__all__ = ['gui_mag_slicer']
_log = logging.getLogger(__name__)
ui_location = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'mag_slicer.ui')
UI_MainWindow, QMainWindow = loadUiType(ui_location)
class Main(QMainWindow, UI_MainWindow):
def __init__(self):
super().__init__()
self.setupUi(self)
self.connect(self.checkBoxLog, QtCore.SIGNAL('clicked()'),
self.update_slice)
self.connect(self.checkBoxScale, QtCore.SIGNAL('clicked()'),
self.update_slice)
self.connect(self.spinBoxSlice, QtCore.SIGNAL('valueChanged(int)'),
self.update_slice)
self.connect(self.comboBoxSlice, QtCore.SIGNAL('currentIndexChanged(int)'),
self.update_phase)
self.connect(self.spinBoxGain, QtCore.SIGNAL('valueChanged(double)'),
self.update_phase)
self.connect(self.checkBoxAuto, QtCore.SIGNAL('toggled(bool)'),
self.update_phase)
self.connect(self.checkBoxSmooth, QtCore.SIGNAL('toggled(bool)'),
self.update_phase)
self.connect(self.pushButtonLoad, QtCore.SIGNAL('clicked()'),
self.load)
self.is_magdata_loaded = False
self.magdata = None
def addmpl(self):
fig = Figure()
fig.add_subplot(111, aspect='equal')
self.canvasMag = FigureCanvas(fig)
self.layoutMag.addWidget(self.canvasMag)
self.canvasMag.draw()
self.toolbarMag = NavigationToolbar(self.canvasMag, self, coordinates=True)
self.layoutMag.addWidget(self.toolbarMag)
fig = Figure()
fig.add_subplot(111, aspect='equal')
self.canvasPhase = FigureCanvas(fig)
self.layoutPhase.addWidget(self.canvasPhase)
self.canvasPhase.draw()
self.toolbarPhase = NavigationToolbar(self.canvasPhase, self, coordinates=True)
self.layoutPhase.addWidget(self.toolbarPhase)
fig = Figure()
fig.add_subplot(111, aspect='equal')
self.canvasHolo = FigureCanvas(fig)
self.layoutHolo.addWidget(self.canvasHolo)
self.canvasHolo.draw()
self.toolbarHolo = NavigationToolbar(self.canvasHolo, self, coordinates=True)
self.layoutHolo.addWidget(self.toolbarHolo)
def update_phase(self):
if self.is_magdata_loaded:
mode_ind = self.comboBoxSlice.currentIndex()
if mode_ind == 0:
self.mode = 'z'
length = self.magdata.dim[0] - 1
elif mode_ind == 1:
self.mode = 'y'
length = self.magdata.dim[1] - 1
else:
self.mode = 'x'
length = self.magdata.dim[2] - 1
if self.checkBoxAuto.isChecked():
gain = 'auto'
else:
gain = self.spinBoxGain.value()
self.projector = SimpleProjector(self.magdata.dim, axis=self.mode)
self.spinBoxSlice.setMaximum(length)
self.scrollBarSlice.setMaximum(length)
self.spinBoxSlice.setValue(int(length / 2.))
self.update_slice()
kernel = Kernel(self.magdata.a, self.projector.dim_uv)
self.phasemapper = PhaseMapperRDFC(kernel)
self.phasemap = self.phasemapper(self.projector(self.magdata))
self.canvasPhase.figure.axes[0].clear()
self.phasemap.plot_phase(axis=self.canvasPhase.figure.axes[0], cbar=False)
if self.checkBoxSmooth.isChecked():
interpolation = 'bilinear'
else:
interpolation = 'none'
self.canvasHolo.figure.axes[0].clear()
self.phasemap.plot_holo(axis=self.canvasHolo.figure.axes[0], gain=gain,
interpolation=interpolation)
self.canvasPhase.draw()
self.canvasHolo.draw()
def update_slice(self):
if self.is_magdata_loaded:
self.canvasMag.figure.axes[0].clear()
self.magdata.plot_quiver(axis=self.canvasMag.figure.axes[0], proj_axis=self.mode,
ax_slice=self.spinBoxSlice.value(),
log=self.checkBoxLog.isChecked(),
scaled=self.checkBoxScale.isChecked())
self.canvasMag.draw()
def load(self):
try:
mag_file = QtGui.QFileDialog.getOpenFileName(self, 'Open Data File', '',
'HDF5 files (*.hdf5)')
except ValueError:
return # Abort if no conf_path is selected!
import hyperspy.api as hs
print(hs.load(mag_file))
self.magdata = load_vectordata(mag_file)
if not self.is_magdata_loaded:
self.addmpl()
self.is_magdata_loaded = True
self.comboBoxSlice.setCurrentIndex(0)
self.update_phase()
def gui_mag_slicer():
"""Call the GUI for viewing magnetic distributions."""
_log.debug('Calling gui_mag_slicer')
app = QtGui.QApplication(sys.argv)
main = Main()
main.show()
app.exec()
return main.magdata
pyramid/utils/mag_slicer.ui deleted 100644 → 0
View file @ f6abe2ad
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>MainWindow</class>
<widget class="QMainWindow" name="MainWindow">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>1222</width>
<height>480</height>
</rect>
</property>
<property name="windowTitle">
<string>MagSlicer</string>
</property>
<widget class="QWidget" name="centralwidget">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Expanding">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="windowTitle">
<string>Mag Slicer</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<item>
<layout class="QHBoxLayout" name="horizontalLayout">
<item>
<widget class="QPushButton" name="pushButtonLoad">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Laden</string>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="checkBoxScale">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Scaled</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
<property name="tristate">
<bool>false</bool>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="checkBoxLog">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Log</string>
</property>
</widget>
</item>
<item>
<widget class="QComboBox" name="comboBoxSlice">
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<item>
<property name="text">
<string>xy-plane</string>
</property>
</item>
<item>
<property name="text">
<string>xz-plane</string>
</property>
</item>
<item>
<property name="text">
<string>zy-plane</string>
</property>
</item>
</widget>
</item>
<item>
<widget class="QLabel" name="labelSlice">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Slice:</string>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
</widget>
</item>
<item>
<widget class="QSpinBox" name="spinBoxSlice">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
<property name="maximum">
<number>0</number>
</property>
<property name="value">
<number>0</number>
</property>
</widget>
</item>
<item>
<widget class="QScrollBar" name="scrollBarSlice">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="minimumSize">
<size>
<width>400</width>
<height>0</height>
</size>
</property>
<property name="maximum">
<number>0</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="checkBoxSmooth">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Smooth</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="checkBoxAuto">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Auto</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
</widget>
</item>
<item>
<widget class="QLabel" name="labelGain">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Gain:</string>
</property>
</widget>
</item>
<item>
<widget class="QDoubleSpinBox" name="spinBoxGain">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximum">
<double>1000000.000000000000000</double>
</property>
<property name="singleStep">
<double>0.100000000000000</double>
</property>
<property name="value">
<double>1.000000000000000</double>
</property>
</widget>
</item>
</layout>
</item>
<item>
<widget class="QWidget" name="layoutPlots" native="true">
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Expanding">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<layout class="QHBoxLayout" name="horizontalLayoutPlots">
<item>
<layout class="QVBoxLayout" name="layoutMag"/>
</item>
<item>
<layout class="QVBoxLayout" name="layoutPhase"/>
</item>
<item>
<layout class="QVBoxLayout" name="layoutHolo"/>
</item>
</layout>
</widget>
</item>
</layout>
</widget>
</widget>
<resources/>
<connections>
<connection>
<sender>spinBoxSlice</sender>
<signal>valueChanged(int)</signal>
<receiver>scrollBarSlice</receiver>
<slot>setValue(int)</slot>
<hints>
<hint type="sourcelabel">
<x>322</x>
<y>19</y>
</hint>
<hint type="destinationlabel">
<x>386</x>
<y>20</y>
</hint>
</hints>
</connection>
<connection>
<sender>scrollBarSlice</sender>
<signal>sliderMoved(int)</signal>
<receiver>spinBoxSlice</receiver>
<slot>setValue(int)</slot>
<hints>
<hint type="sourcelabel">
<x>431</x>
<y>24</y>
</hint>
<hint type="destinationlabel">
<x>321</x>
<y>20</y>
</hint>
</hints>
</connection>
<connection>
<sender>checkBoxAuto</sender>
<signal>toggled(bool)</signal>
<receiver>spinBoxGain</receiver>
<slot>setDisabled(bool)</slot>
<hints>
<hint type="sourcelabel">
<x>1065</x>
<y>22</y>
</hint>
<hint type="destinationlabel">
<x>1166</x>
<y>21</y>
</hint>
</hints>
</connection>
</connections>
</ui>
pyramid/utils/phasemap_creator.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
# Form implementation generated from reading ui file 'phasemap_creator.ui'
#
# Created: Thu Sep 24 11:42:11 2015
# by: PyQt4 UI code generator 4.9.6
#
# WARNING! All changes made in this file will be lost!
"""GUI for setting up PhasMaps from existing data in different formats."""
import logging
import os
import sys
from matplotlib.figure import Figure
from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
from PIL import Image
import numpy as np
#import hyperspy.api as hs # TODO: Necessary?
import pyramid as pr
try:
from PyQt5 import QtGui, QtCore
from PyQt5.uic import loadUiType
except ImportError:
from PyQt4 import QtGui, QtCore
from PyQt4.uic import loadUiType
__all__ = ['gui_phasemap_creator']
_log = logging.getLogger(__name__)
ui_location = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'phasemap_creator.ui')
UI_MainWindow, QMainWindow = loadUiType(ui_location)
class Main(QMainWindow, UI_MainWindow):
def __init__(self):
super().__init__()
self.setupUi(self)
self.connect(self.pushButton_phase, QtCore.SIGNAL('clicked()'),
self.load_phase)
self.connect(self.pushButton_mask, QtCore.SIGNAL('clicked()'),
self.load_mask)
self.connect(self.pushButton_conf, QtCore.SIGNAL('clicked()'),
self.load_conf)
self.connect(self.pushButton_export, QtCore.SIGNAL('clicked()'),
self.export)
self.connect(self.horizontalScrollBar, QtCore.SIGNAL('valueChanged(int)'),
self.doubleSpinBox_thres.setValue)
self.connect(self.doubleSpinBox_thres, QtCore.SIGNAL('valueChanged(double)'),
self.horizontalScrollBar.setValue)
self.connect(self.checkBox_mask, QtCore.SIGNAL('clicked()'),
self.update_phasemap)
self.connect(self.checkBox_conf, QtCore.SIGNAL('clicked()'),
self.update_phasemap)
self.connect(self.doubleSpinBox_a, QtCore.SIGNAL('editingFinished()'),
self.update_phasemap)
self.connect(self.doubleSpinBox_thres, QtCore.SIGNAL('valueChanged(double)'),
self.update_mask)
self.phase_loaded = False
self.mask_loaded = False
self.dir = ''
self.phasemap = None
def addmpl(self):
fig = Figure()
fig.add_subplot(111, aspect='equal')
self.canvas = FigureCanvas(fig)
self.mplLayout.addWidget(self.canvas)
self.canvas.draw()
self.toolbar = NavigationToolbar(self.canvas, self, coordinates=True)
self.mplLayout.addWidget(self.toolbar)
def update_phasemap(self):
if self.phase_loaded:
self.phasemap.a = self.doubleSpinBox_a.value()
show_mask = self.checkBox_mask.isChecked()
show_conf = self.checkBox_conf.isChecked()
self.canvas.figure.axes[0].clear()
self.canvas.figure.axes[0].hold(True)
self.phasemap.plot_phase('PhaseMap', axis=self.canvas.figure.axes[0],
show_mask=show_mask, show_conf=show_conf, cbar=False)
self.canvas.draw()
def update_mask(self):
if self.mask_loaded:
threshold = self.doubleSpinBox_thres.value()
mask_img = Image.fromarray(self.raw_mask)
mask = np.asarray(mask_img.resize(list(reversed(self.phasemap.dim_uv))))
self.phasemap.mask = np.where(mask >= threshold, True, False)
self.update_phasemap()
def load_phase(self):
try:
self.phase_path = QtGui.QFileDialog.getOpenFileName(self, 'Load Phase', self.dir)
self.phasemap = pr.file_io.io_phasemap._load(self.phase_path, as_phasemap=True)
except ValueError:
return # Abort if no phase_path is selected!
self.doubleSpinBox_a.setValue(self.phasemap.a)
self.dir = os.path.join(os.path.dirname(self.phase_path))
if not self.phase_loaded:
self.addmpl()
self.pushButton_mask.setEnabled(True)
self.pushButton_conf.setEnabled(True)
self.pushButton_export.setEnabled(True)
self.phase_loaded = True
self.horizontalScrollBar.setMinimum(0)
self.horizontalScrollBar.setMaximum(0)
self.horizontalScrollBar.setEnabled(False)
self.doubleSpinBox_thres.setMinimum(0)
self.doubleSpinBox_thres.setMaximum(0)
self.doubleSpinBox_thres.setValue(0)
self.doubleSpinBox_thres.setEnabled(False)
self.mask_loaded = False
self.update_phasemap()
def load_mask(self):
try:
mask_path = QtGui.QFileDialog.getOpenFileName(self, 'Load Mask', self.dir)
self.raw_mask = pr.file_io.io_phasemap._load(mask_path)
except ValueError:
return # Abort if no mask_path is selected!
mask_min = self.raw_mask.min()
mask_max = self.raw_mask.max()
self.horizontalScrollBar.setEnabled(True)
self.horizontalScrollBar.setMinimum(mask_min)
self.horizontalScrollBar.setMaximum(mask_max)
self.horizontalScrollBar.setSingleStep((mask_max - mask_min) / 255.)
self.horizontalScrollBar.setValue((mask_max - mask_min) / 2.)
self.doubleSpinBox_thres.setEnabled(True)
self.doubleSpinBox_thres.setMinimum(mask_min)
self.doubleSpinBox_thres.setMaximum(mask_max)
self.doubleSpinBox_thres.setSingleStep((mask_max - mask_min) / 255.)
self.doubleSpinBox_thres.setValue((mask_max - mask_min) / 2.)
self.mask_loaded = True
self.update_mask()
def load_conf(self):
try:
conf_path = QtGui.QFileDialog.getOpenFileName(self, 'Load Confidence', self.dir)
confidence = pr.file_io.io_phasemap._load(conf_path)
except ValueError:
return # Abort if no conf_path is selected!
confidence = confidence.astype(float) / confidence.max() + 1e-30
self.phasemap.confidence = confidence
self.update_phasemap()
def export(self):
try:
export_name = os.path.splitext(os.path.basename(self.phase_path))[0]
export_default = os.path.join(self.dir, 'phasemap_gui_{}.hdf5'.format(export_name))
export_path = QtGui.QFileDialog.getSaveFileName(self, 'Export PhaseMap',
export_default, 'HDF5 (*.hdf5)')
self.phasemap.to_signal().save(export_path, overwrite=True)
except (ValueError, AttributeError):
return # Abort if no export_path is selected or self.phasemap doesn't exist yet!
def gui_phasemap_creator():
"""Call the GUI for phasemap creation."""
_log.debug('Calling gui_phasemap_creator')
app = QtGui.QApplication(sys.argv)
main = Main()
main.show()
app.exec()
return main.phasemap
pyramid/utils/phasemap_creator.ui deleted 100644 → 0
View file @ f6abe2ad
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>MainWindow</class>
<widget class="QMainWindow" name="MainWindow">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>726</width>
<height>632</height>
</rect>
</property>
<property name="windowTitle">
<string>MainWindow</string>
</property>
<widget class="QWidget" name="centralwidget">
<property name="windowTitle">
<string>Form</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<item>
<widget class="QWidget" name="mplwidget" native="true">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Expanding">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<layout class="QVBoxLayout" name="verticalLayout_2">
<item>
<layout class="QVBoxLayout" name="mplLayout"/>
</item>
</layout>
</widget>
</item>
<item>
<layout class="QHBoxLayout" name="horizontalLayout">
<item>
<widget class="QPushButton" name="pushButton_phase">
<property name="text">
<string>Load Phase</string>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="pushButton_mask">
<property name="text">
<string>Load Mask</string>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="pushButton_conf">
<property name="text">
<string>Load Confidence</string>
</property>
<property name="checkable">
<bool>false</bool>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="pushButton_export">
<property name="text">
<string>Export Phasemap</string>
</property>
</widget>
</item>
</layout>
</item>
<item>
<layout class="QHBoxLayout" name="horizontalLayout_2">
<item>
<widget class="QLabel" name="label_2">
<property name="text">
<string>Grid spacing [nm]:</string>
</property>
</widget>
</item>
<item>
<widget class="QDoubleSpinBox" name="doubleSpinBox_a">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="minimumSize">
<size>
<width>60</width>
<height>0</height>
</size>
</property>
<property name="baseSize">
<size>
<width>0</width>
<height>0</height>
</size>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
<property name="maximum">
<double>1000.000000000000000</double>
</property>
<property name="value">
<double>1.000000000000000</double>
</property>
</widget>
</item>
<item>
<widget class="QLabel" name="label">
<property name="sizePolicy">
<sizepolicy hsizetype="Fixed" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Mask Threshold:</string>
</property>
</widget>
</item>
<item>
<widget class="QDoubleSpinBox" name="doubleSpinBox_thres">
<property name="sizePolicy">
<sizepolicy hsizetype="Minimum" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="minimumSize">
<size>
<width>60</width>
<height>0</height>
</size>
</property>
<property name="maximumSize">
<size>
<width>16777215</width>
<height>16777215</height>
</size>
</property>
<property name="baseSize">
<size>
<width>0</width>
<height>0</height>
</size>
</property>
<property name="alignment">
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
<property name="prefix">
<string/>
</property>
<property name="decimals">
<number>2</number>
</property>
<property name="maximum">
<double>0.000000000000000</double>
</property>
</widget>
</item>
<item>
<widget class="QScrollBar" name="horizontalScrollBar">
<property name="sizePolicy">
<sizepolicy hsizetype="MinimumExpanding" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximum">
<number>0</number>
</property>
<property name="singleStep">
<number>1</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="checkBox_mask">
<property name="sizePolicy">
<sizepolicy hsizetype="Fixed" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>show mask</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="checkBox_conf">
<property name="sizePolicy">
<sizepolicy hsizetype="Fixed" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>show confidence</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
</widget>
</item>
</layout>
</item>
</layout>
</widget>
</widget>
<resources/>
<connections/>
</ui>
pyramid/utils/pm.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Convenience function for phase mapping magnetic or charge distributions."""
import logging
from ..kernel import Kernel, KernelCharge
from ..phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperCharge
from ..projector import RotTiltProjector, XTiltProjector, YTiltProjector, SimpleProjector
__all__ = ['pm']
_log = logging.getLogger(__name__)
def pm(fielddata, mode='z', b_0=1, electrode_vec=(1E6, 1E6), mapper='RDFC', **kwargs):
"""Convenience function for fast magnetic phase mapping.
Parameters
----------
fielddata : :class:`~.VectorData`, or `~.ScalarData`
A :class:`~.VectorData` or `~.ScalarData` object, from which the projected phase map should be calculated.
mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
is used for the projection. Default is a simple projection along the `z`-direction.
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
electrode_vec : tuple of float (N=2)
The norm vector of the counter electrode, (elec_a,elec_b), and the distance to the origin is
the norm of (elec_a,elec_b). The default value is (1E6, 1E6).
mapper : :class: '~. PhaseMap'
A :class: '~. PhaseMap' object, which maps a fielddata into a phase map. The default is 'RDFC'.
**kwargs : additional arguments
Additional arguments like `dim_uv`, 'tilt' or 'rotation', which are passed to the
projector-constructor, defined by the `mode`.
Returns
-------
phasemap : :class:`~pyramid.phasemap.PhaseMap`
The calculated phase map as a :class:`~.PhaseMap` object.
"""
_log.debug('Calling pm')
# In case of FDFC:
padding = kwargs.pop('padding', 0)
# Determine projection mode:
if mode == 'rot-tilt':
projector = RotTiltProjector(fielddata.dim, **kwargs)
elif mode == 'x-tilt':
projector = XTiltProjector(fielddata.dim, **kwargs)
elif mode == 'y-tilt':
projector = YTiltProjector(fielddata.dim, **kwargs)
elif mode in ['x', 'y', 'z']:
projector = SimpleProjector(fielddata.dim, axis=mode, **kwargs)
else:
raise ValueError("Invalid mode (use 'x', 'y', 'z', 'x-tilt', 'y-tilt' or 'rot-tilt')")
# Project:
field_proj = projector(fielddata)
# Set up phasemapper and map phase:
if mapper == 'RDFC':
phasemapper = PhaseMapperRDFC(Kernel(fielddata.a, projector.dim_uv, b_0=b_0))
elif mapper == 'FDFC':
phasemapper = PhaseMapperFDFC(fielddata.a, projector.dim_uv, b_0=b_0, padding=padding)
# Set up phasemapper and map phase:
if mapper == 'Charge':
phasemapper = PhaseMapperCharge(KernelCharge(fielddata.a, projector.dim_uv, electrode_vec=electrode_vec))
else:
raise ValueError("Invalid mapper (use 'RDFC', 'FDFC' or 'Charge'")
phasemap = phasemapper(field_proj)
# Get mask from fielddata:
phasemap.mask = field_proj.get_mask()[0, ...]
# Return phase:
return phasemap
def pm2(magdata, mode='z', b_0=1, mapper='RDFC', **kwargs):
"""Convenience function for fast magnetic phase mapping.
Parameters
----------
magdata : :class:`~.VectorData`
A :class:`~.VectorData` object, from which the projected phase map should be calculated.
mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
is used for the projection. Default is a simple projection along the `z`-direction.
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
**kwargs : additional arguments
Additional arguments like `dim_uv`, 'tilt' or 'rotation', which are passed to the
projector-constructor, defined by the `mode`.
Returns
-------
phasemap : :class:`~pyramid.phasemap.PhaseMap`
The calculated phase map as a :class:`~.PhaseMap` object.
"""
_log.debug('Calling pm2')
# In case of FDFC:
padding = kwargs.pop('padding', 0)
# Determine projection mode:
if mode == 'rot-tilt':
projector = RotTiltProjector(magdata.dim, **kwargs)
elif mode == 'x-tilt':
projector = XTiltProjector(magdata.dim, **kwargs)
elif mode == 'y-tilt':
projector = YTiltProjector(magdata.dim, **kwargs)
elif mode in ['x', 'y', 'z']:
projector = SimpleProjector(magdata.dim, axis=mode, **kwargs)
else:
raise ValueError("Invalid mode (use 'x', 'y', 'z', 'x-tilt', 'y-tilt' or 'rot-tilt')")
# Project:
mag_proj = projector(magdata)
# Set up phasemapper and map phase:
if mapper == 'RDFC':
phasemapper = PhaseMapperRDFC(Kernel(magdata.a, projector.dim_uv, b_0=b_0))
elif mapper == 'FDFC':
phasemapper = PhaseMapperFDFC(magdata.a, projector.dim_uv, b_0=b_0, padding=padding)
else:
raise ValueError("Invalid mapper (use 'RDFC' or 'FDFC'")
phasemap = phasemapper(mag_proj)
# Get mask from magdata:
phasemap.mask = mag_proj.get_mask()[0, ...]
# Return phase:
return phasemap
def pm3(elecdata, mode='z', electrode_vec=(1E6, 1E6), mapper='Charge', **kwargs):
"""Convenience function for fast electric phase mapping.
Parameters
----------
elecdata : :class:`~.ScalarData`
A :class:`~.ScalarData` object, from which the projected phase map should be calculated.
mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
is used for the projection. Default is a simple projection along the `z`-direction.
electrode_vec : tuple of float (N=2)
The norm vector of the counter electrode, (elec_a,elec_b), and the distance to the origin is
the norm of (elec_a,elec_b).
**kwargs : additional arguments
Additional arguments like `dim_uv`, 'tilt' or 'rotation', which are passed to the
projector-constructor, defined by the `mode`.
Returns
-------
phasemap : :class:`~pyramid.phasemap.PhaseMap`
The calculated phase map as a :class:`~.PhaseMap` object.
"""
_log.debug('Calling pm3')
# Determine projection mode:
if mode == 'rot-tilt':
projector = RotTiltProjector(elecdata.dim, **kwargs)
elif mode == 'x-tilt':
projector = XTiltProjector(elecdata.dim, **kwargs)
elif mode == 'y-tilt':
projector = YTiltProjector(elecdata.dim, **kwargs)
elif mode in ['x', 'y', 'z']:
projector = SimpleProjector(elecdata.dim, axis=mode, **kwargs)
else:
raise ValueError("Invalid mode (use 'x', 'y', 'z', 'x-tilt', 'y-tilt' or 'rot-tilt')")
# Project:
charge_proj = projector(elecdata)
# Set up phasemapper and map phase:
if mapper == 'Charge':
phasemapper = PhaseMapperCharge(KernelCharge(elecdata.a, projector.dim_uv, electrode_vec=electrode_vec))
else:
raise ValueError("Invalid mapper (use 'Charge'")
phasemap = phasemapper(charge_proj)
# Get mask from elecdata:
phasemap.mask = charge_proj.get_mask()[0, ...]
# Return phase:
return phasemap
\ No newline at end of file
pyramid/utils/reconstruction_2d_from_phasemap.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Reconstruct a magnetization distributions from a single phase map."""
import logging
import numpy as np
from .. import reconstruction
from ..dataset import DataSet
from ..projector import SimpleProjector
from ..regularisator import FirstOrderRegularisator, NoneRegularisator
from ..kernel import KernelCharge
from ..phasemapper import PhaseMapperCharge
from ..forwardmodel import ForwardModel
from ..costfunction import Costfunction
from .pm import pm
__all__ = ['reconstruction_2d_from_phasemap']
_log = logging.getLogger(__name__)
def reconstruction_2d_from_phasemap(phasemap, b_0=1, lam=1E-3, max_iter=100, ramp_order=1,
plot_results=False, ar_dens=None, verbose=True):
"""Convenience function for reconstructing a projected distribution from a single phasemap.
Parameters
----------
phasemap: :class:`~PhaseMap`
The phasemap which is used for the reconstruction.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
lam : float
Regularisation parameter determining the weighting between measurements and regularisation.
max_iter : int, optional
The maximum number of iterations for the opimization.
ramp_order : int or None (default)
Polynomial order of the additional phase ramp which will be added to the phase maps.
All ramp parameters have to be at the end of the input vector and are split automatically.
Default is None (no ramps are added).
plot_results: boolean, optional
If True, the results are plotted after reconstruction.
ar_dens: int, optional
Number defining the arrow density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
verbose: bool, optional
If set to True, information like a progressbar is displayed during reconstruction.
The default is False.
Returns
-------
magdata_rec, cost: :class:`~.VectorData`, :class:`~.Costfunction`
The reconstructed magnetisation distribution and the used costfunction.
"""
_log.debug('Calling reconstruction_2d_from_phasemap')
# Construct DataSet, Regularisator, ForwardModel and Costfunction:
dim = (1,) + phasemap.dim_uv
data = DataSet(phasemap.a, dim, b_0)
data.append(phasemap, SimpleProjector(dim))
data.set_3d_mask()
fwd_model = ForwardModel(data, ramp_order)
reg = FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
cost = Costfunction(fwd_model, reg)
# Reconstruct:
magdata_rec = reconstruction.optimize_linear(cost, max_iter=max_iter, verbose=verbose)
param_cache = cost.fwd_model.ramp.param_cache
if ramp_order is None:
offset, ramp = 0, (0, 0)
elif ramp_order >= 1:
offset, ramp = param_cache[0][0], (param_cache[1][0], param_cache[2][0])
elif ramp_order == 0:
offset, ramp = param_cache[0][0], (0, 0)
else:
raise ValueError('ramp_order has to be a positive integer or None!')
# Plot stuff:
if plot_results:
if ar_dens is None:
ar_dens = np.max([1, np.max(dim) // 64])
magdata_rec.plot_quiver_field(note='Reconstructed Distribution',
ar_dens=ar_dens, figsize=(16, 16))
phasemap_rec = pm(magdata_rec)
gain = 4 * 2 * np.pi / (np.abs(phasemap_rec.phase).max() + 1E-30)
gain = round(gain, -int(np.floor(np.log10(abs(gain)))))
vmin = phasemap_rec.phase.min()
vmax = phasemap_rec.phase.max()
phasemap.plot_combined(note='Input Phase', gain=gain)
phasemap -= fwd_model.ramp(index=0)
phasemap.plot_combined(note='Input Phase (ramp corrected)', gain=gain, vmin=vmin, vmax=vmax)
title = 'Reconstructed Phase'
if ramp_order is not None:
if ramp_order >= 0:
print('offset:', offset)
# title += ', fitted Offset: {:.2g} [rad]'.format(offset)
if ramp_order >= 1:
print('ramp:', ramp)
# title += ', (Fitted Ramp: (u:{:.2g}, v:{:.2g}) [rad/nm]'.format(*ramp)
phasemap_rec.plot_combined(note=title, gain=gain, vmin=vmin, vmax=vmax)
diff = (phasemap_rec - phasemap)
diff_name = 'Difference (RMS: {:.2g} rad)'.format(np.sqrt(np.mean(diff.phase) ** 2))
diff.plot_phase_with_hist(note=diff_name, sigma_clip=3)
if ramp_order is not None:
ramp = fwd_model.ramp(0)
ramp.plot_phase(note='Fitted Ramp')
# Return reconstructed magnetisation distribution and cost function:
return magdata_rec, cost
def reconstruction_2d_charge_from_phasemap(phasemap, lam=1E-3, max_iter=100, ramp_order=1,
electrode_vec=(1E6, 1E6), v_acc=300E3,
plot_results=False, verbose=True):
"""Convenience function for reconstructing a projected distribution from a single phasemap.
Parameters
----------
phasemap: :class:`~PhaseMap`
The phasemap which is used for the reconstruction.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
lam : float
Regularisation parameter determining the weighting between measurements and regularisation.
max_iter : int, optional
The maximum number of iterations for the opimization.
ramp_order : int or None (default)
Polynomial order of the additional phase ramp which will be added to the phase maps.
All ramp parameters have to be at the end of the input vector and are split automatically.
Default is None (no ramps are added).
plot_results: boolean, optional
If True, the results are plotted after reconstruction.
ar_dens: int, optional
Number defining the arrow density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
verbose: bool, optional
If set to True, information like a progressbar is displayed during reconstruction.
The default is False.
Returns
-------
magdata_rec, cost: :class:`~.VectorData`, :class:`~.Costfunction`
The reconstructed magnetisation distribution and the used costfunction.
"""
_log.debug('Calling reconstruction_2d_from_phasemap')
# Construct DataSet, Regularisator, ForwardModel and Costfunction:
dim = (1,) + phasemap.dim_uv
data = DataSet(phasemap.a, dim, b_0=1)
kernel = KernelCharge(phasemap.a, phasemap.dim_uv, electrode_vec=electrode_vec, v_acc=v_acc)
data.append(phasemap, SimpleProjector(dim), PhaseMapperCharge(kernel))
data.set_3d_mask()
# TODO: Rework classes below (ForwardModel, Costfunction)!
fwd_model = ForwardModel(data, ramp_order)
reg = NoneRegularisator()#FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
cost = Costfunction(fwd_model, reg)
# Reconstruct:
magdata_rec = reconstruction.optimize_linear(cost, max_iter=max_iter, verbose=verbose)
param_cache = cost.fwd_model.ramp.param_cache
if ramp_order is None:
offset, ramp = 0, (0, 0)
elif ramp_order >= 1:
offset, ramp = param_cache[0][0], (param_cache[1][0], param_cache[2][0])
elif ramp_order == 0:
offset, ramp = param_cache[0][0], (0, 0)
else:
raise ValueError('ramp_order has to be a positive integer or None!')
# Plot stuff:
if plot_results:
if ar_dens is None:
ar_dens = np.max([1, np.max(dim) // 64])
magdata_rec.plot_quiver_field(note='Reconstructed Distribution',
ar_dens=ar_dens, figsize=(16, 16))
phasemap_rec = pm(magdata_rec)
gain = 4 * 2 * np.pi / (np.abs(phasemap_rec.phase).max() + 1E-30)
gain = round(gain, -int(np.floor(np.log10(abs(gain)))))
vmin = phasemap_rec.phase.min()
vmax = phasemap_rec.phase.max()
phasemap.plot_combined(note='Input Phase', gain=gain)
phasemap -= fwd_model.ramp(index=0)
phasemap.plot_combined(note='Input Phase (ramp corrected)', gain=gain, vmin=vmin, vmax=vmax)
title = 'Reconstructed Phase'
if ramp_order is not None:
if ramp_order >= 0:
print('offset:', offset)
# title += ', fitted Offset: {:.2g} [rad]'.format(offset)
if ramp_order >= 1:
print('ramp:', ramp)
# title += ', (Fitted Ramp: (u:{:.2g}, v:{:.2g}) [rad/nm]'.format(*ramp)
phasemap_rec.plot_combined(note=title, gain=gain, vmin=vmin, vmax=vmax)
diff = (phasemap_rec - phasemap)
diff_name = 'Difference (RMS: {:.2g} rad)'.format(np.sqrt(np.mean(diff.phase) ** 2))
diff.plot_phase_with_hist(note=diff_name, sigma_clip=3)
if ramp_order is not None:
ramp = fwd_model.ramp(0)
ramp.plot_phase(note='Fitted Ramp')
# Return reconstructed magnetisation distribution and cost function:
return magdata_rec, cost
pyramid/utils/reconstruction_3d_from_magdata.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Reconstruct a magnetization distributions from phase maps created from it."""
import logging
import numpy as np
import multiprocessing as mp
from .. import reconstruction
from ..dataset import DataSet
from ..projector import XTiltProjector, YTiltProjector
from ..ramp import Ramp
from ..regularisator import FirstOrderRegularisator
from ..forwardmodel import ForwardModel, DistributedForwardModel
from ..costfunction import Costfunction
from ..phasemapper import PhaseMapperRDFC
from ..kernel import Kernel
__all__ = ['reconstruction_3d_from_magdata']
_log = logging.getLogger(__name__)
def reconstruction_3d_from_magdata(magdata, b_0=1, lam=1E-3, max_iter=100, ramp_order=1,
angles=np.linspace(-90, 90, num=19), dim_uv=None,
axes=(True, True), noise=0, offset_max=0, ramp_max=0,
use_internal_mask=True, plot_results=False, plot_input=False,
ar_dens=None, multicore=False, verbose=True):
"""Convenience function for reconstructing a projected distribution from a single phasemap.
Parameters
----------
magdata: :class:`~.VectorData`
The magnetisation distribution which should be used for the reconstruction.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
lam : float
Regularisation parameter determining the weighting between measurements and regularisation.
max_iter : int, optional
The maximum number of iterations for the opimization.
ramp_order : int or None (default)
Polynomial order of the additional phase ramp which will be added to the phase maps.
All ramp parameters have to be at the end of the input vector and are split automatically.
Default is None (no ramps are added).
angles: :class:`~numpy.ndarray` (N=1), optional
Numpy array determining the angles which should be used for the projectors in x- and
y-direction. This implicitly sets the number of images per rotation axis. Defaults to a
range from -90° to 90° degrees, in 10° steps.
dim_uv: int or None (default)
Determines if the phasemaps should be padded to a certain size while calculating.
axes: tuple of booleans (N=2), optional
Determines if both tilt axes should be calculated. The order is (x, y), both are True by
default.
noise: float, optional
If this is not zero, random gaussian noise with this as a maximum value will be applied
to all calculated phasemaps. The default is 0. The unit is radians.
offset_max: float, optional
if this is not zero, a random offset with this as a maximum value will be applied to all
calculated phasemaps. The default is 0.
ramp_max: float, optional
if this is not zero, a random linear ramp with this as a maximum value will be applied
to both axes of all calculated phasemaps. The default is 0.
use_internal_mask: boolean, optional
If True, the mask from the input magnetization distribution is taken for the
reconstruction. If False, the mask is calculated via logic backprojection from the 2D-masks
of the input phasemaps.
plot_results: boolean, optional
If True, the results are plotted after reconstruction.
plot_input:
If True, the input phasemaps are plotted after reconstruction.
ar_dens: int, optional
Number defining the arrow density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
multicore: boolean, optional
Determines if multiprocessing should be used. Default is True. Phasemap calculations
will be divided onto the separate cores.
verbose: bool, optional
If set to True, information like a progressbar is displayed during reconstruction.
The default is False.
Returns
-------
magdata_rec, cost: :class:`~.VectorData`, :class:`~.Costfunction`
The reconstructed magnetisation distribution and the used costfunction.
"""
_log.debug('Calling reconstruction_3d_from_magdata')
# Construct DataSet:
dim = magdata.dim
if ar_dens is None:
ar_dens = np.max([1, np.max(dim) // 128])
data = DataSet(magdata.a, magdata.dim, b_0)
# Construct projectors:
projectors = []
# Construct data set and regularisator:
for angle in angles:
angle_rad = angle * np.pi / 180
if axes[0]:
projectors.append(XTiltProjector(magdata.dim, angle_rad, dim_uv))
if axes[1]:
projectors.append(YTiltProjector(magdata.dim, angle_rad, dim_uv))
# Add pairs of projectors and according phasemaps to the DataSet:
for projector in projectors:
mag_proj = projector(magdata)
phasemap = PhaseMapperRDFC(Kernel(magdata.a, projector.dim_uv, b_0))(mag_proj)
phasemap.mask = mag_proj.get_mask()[0, ...]
data.append(phasemap, projector)
# Add offset and ramp if necessary:
for i, phasemap in enumerate(data.phasemaps):
offset = np.random.uniform(-offset_max, offset_max)
ramp_u = np.random.uniform(-ramp_max, ramp_max)
ramp_v = np.random.uniform(-ramp_max, ramp_max)
phasemap += Ramp.create_ramp(phasemap.a, phasemap.dim_uv, (offset, ramp_u, ramp_v))
data.phasemaps[i] = phasemap
# Add noise if necessary:
if noise != 0: # TODO: write function to add noise after APERTURE!! (ask Florian again)
for i, phasemap in enumerate(data.phasemaps):
phasemap.phase += np.random.normal(0, noise, phasemap.dim_uv)
data.phasemaps[i] = phasemap
# Construct mask:
if use_internal_mask:
data.mask = magdata.get_mask() # Use perfect mask from magdata!
else:
data.set_3d_mask() # Construct mask from 2D phase masks!
# Construct regularisator, forward model and costfunction:
if multicore:
mp.freeze_support()
fwd_model = DistributedForwardModel(data, ramp_order=ramp_order, nprocs=mp.cpu_count())
else:
fwd_model = ForwardModel(data, ramp_order=ramp_order)
reg = FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
cost = Costfunction(fwd_model, reg)
# Reconstruct and save:
magdata_rec = reconstruction.optimize_linear(cost, max_iter=max_iter, verbose=verbose)
# Finalize ForwardModel (returns workers if multicore):
fwd_model.finalize()
# Plot input:
if plot_input:
data.plot_phasemaps()
# Plot results:
if plot_results:
data.plot_mask()
magdata.plot_quiver3d('Original Distribution', ar_dens=ar_dens)
magdata_rec.plot_quiver3d('Reconstructed Distribution (angle)', ar_dens=ar_dens)
magdata_rec.plot_quiver3d('Reconstructed Distribution (amplitude)',
ar_dens=ar_dens, coloring='amplitude')
# Return reconstructed magnetisation distribution and cost function:
return magdata_rec, cost
pyramid/version.py deleted 100644 → 0
View file @ f6abe2ad
# -*- coding: utf-8 -*-
""""This file is generated automatically by the Pyramid `setup.py`"""
version = "0.1.0-dev"
hg_revision = "???"
requirements.txt deleted 100644 → 0
View file @ f6abe2ad
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
setup.cfg
View file @ 6e3128b8
# -*- coding: utf-8 -*-
#
# mslib.setup
# ~~~~~~~~~~~
#
# setup.cfg
#
# This file is part of mss.
#
# :copyright: Copyright 2016-2017 Reimar Bauer, Joern Ungermann
# :copyright: Copyright 2016-2017 by the mss team, see AUTHORS.
# :license: APACHE-2.0, see LICENSE for details.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in com#pliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# setup.cfg
# :copyright: Copyright 2020 Jan Caron
# TODO: Check if everything is taylored to pyramid!
# TODO: pip install .[io] is not working because hyperspy depends on trait which has no wheels on PyPI at the moment...
# TODO: See https://github.com/hyperspy/hyperspy/issues/2315 and https://github.com/enthought/traits/issues/357
# TODO: Add hyperspy back as soon (if?) this is resolved... Until then: install hyperspy with conda!
# TODO: mayavi has the same exact problem and tvtk is apparently also a problem...
# CONFIGURATION FOR TESTING:
[aliases]
test=pytest
test = pytest
[coverage:run]
branch = True
source = pyramid
source = src/empyre
omit =
tests/*
[tool:pytest]
addopts = --cov --flake8
flake8-max-line-length = 100
flake8-ignore =
ALL # TODO: PEP8 deactivated by this line (remove at a later point)!
E402 E124 E125
pyramid/__init__.py F401
doc/conf.py ALL
scripts/*.py ALL
[flake8]
max-line-length = 120
ignore =
# module import not at top of file
E402
# closing bracket does not match visual indentation
E124
# continuation line with same indent as next logical line
E125
# missing whitespace around arithmetic operator
E226
# line break before binary operator
W503
# line break after binary operator
W504
# do not use variables named ‘l’, ‘O’, or ‘I’
E741
per-file-ignores =
*/__init__.py: F401, F403, F405, F821
# F401: module imported but unused
# F403: 'from module import *' used; unable to detect undefined names
# F405: Name may be undefined, or defined from star imports: module
# F821: undefined name 'name'
setup.py deleted 100644 → 0
View file @ f6abe2ad
#!python
# coding=utf-8
"""Setup for testing, building, distributing and installing the 'Pyramid'-package"""
import os
import re
import subprocess
import sys
import itertools
#from distutils.command.build import build
#import numpy
from setuptools import setup, find_packages
DISTNAME = 'pyramid'
DESCRIPTION = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions'
MAINTAINER = 'Jan Caron'
MAINTAINER_EMAIL = 'j.caron@fz-juelich.de'
URL = ''
VERSION = '0.1.0.dev0' # TODO: Better way?
PYTHON_VERSION = (2, 7) # TODO: get rid of!!!
DEPENDENCIES = {'numpy': (1, 10)} # TODO: get rid of!!!
LONG_DESCRIPTION = 'long description (TODO!)' # TODO: Long description! put in (Readme?) file!
# TODO: get rid of superfluous functions!
def get_package_version(package):
"""Return the package version of the specified package.
Parameters
----------
package: basestring
Name of the package whic should be checked.
Returns
-------
version: tuple (N=3)
Version number as a tuple.
"""
version = []
for version_attr in ('version', 'VERSION', '__version__'):
if (hasattr(package, version_attr) and
isinstance(getattr(package, version_attr), str)):
version_info = getattr(package, version_attr, '')
for part in re.split('\D+', version_info):
try:
version.append(int(part))
except ValueError:
pass
return tuple(version)
def check_requirements():
"""Checks the requirements of the Pyramid package."""
if sys.version_info < PYTHON_VERSION:
raise SystemExit('You need Python version %d.%d or later.'
% PYTHON_VERSION)
for package_name, min_version in DEPENDENCIES.items():
dep_error = False
try:
package = __import__(package_name)
except ImportError:
dep_error = True
else:
package_version = get_package_version(package)
if min_version > package_version:
dep_error = True
if dep_error:
raise ImportError('You need `%s` version %d.%d or later.'
% ((package_name,) + min_version))
def hg_version(): # TODO: Replace with GIT! Also check build output on GitLab! See numpy setup.py!
"""Get the Mercurial reference identifier.
Returns
-------
hg_ref: basestring
The Mercurial reference identifier.
"""
try:
hg_rev = subprocess.check_output(['hg', 'id', '--id']).strip()
except:
hg_rev = "???"
return hg_rev
def write_version_py(filename='pyramid/version.py'):
"""Write the version.py file.
Parameters
----------
filename: basestring, optional
Write the version and hg_revision into the specified python file.
Defaults to 'pyramid/version.py'.
"""
version_string = '# -*- coding: utf-8 -*-\n' + \
'""""This file is generated automatically by the Pyramid `setup.py`"""\n' + \
'version = "{}"\n'.format(VERSION) + \
'hg_revision = "{}"\n'.format(hg_version())
with open(os.path.join(os.path.dirname(__file__), filename), 'w') as vfile:
vfile.write(version_string)
def get_files(rootdir):
"""Returns a list of .py-files inside rootdir.
Parameters
----------
rootdir: basestring
Root directory in which to search for ``.py``-files.
Returns
-------
filepaths: list
List of filepaths which were found.
"""
filepaths = []
for root, dirs, files in os.walk(rootdir):
for filename in files:
if filename.endswith('.py'):
filepaths.append(os.path.join(root, filename))
return filepaths
# TODO: Outsource stuff to setup.cfg? See https://github.com/pypa/setuptools/pull/862
# TODO: Use requirements.txt? extras_require for optional stuff (hyperspy, plotting)?
# TODO: After split of Pyramid, comment out and see what really is used (is e.g. scipy?)!
install_requires = ['numpy>=1.6', 'tqdm', 'scipy', 'matplotlib', 'Pillow', 'h5py',
'hyperspy', 'jutil', 'cmocean']
# TODO: extend extras_require for plotting and IO:
# TODO: extra: 'pyfftw', 'mayavi' (not easy to install... find a way!)
# TODO: See https://stackoverflow.com/a/28842733 for extras_require...
# TODO: ...replace [dev] with [IO] (hyperspy) and [plotting] (separate plotting library)!
extras_require = {
# TODO: Test all if really needed! don't use nose, if possible (pure pytest)!
'tests': ['pytest', 'pytest-runner', 'pytest-cov', 'pytest-flake8', 'coverage'],
'3Dplot': ['qt==4.8', 'mayavi==4.5']
# TODO: more for mayavi (plotting in general) and hyperspy, etc (see below)...
}
extras_require["all"] = list(itertools.chain(*list(extras_require.values())))
print('\n-------------------------------------------------------------------------------')
# print('checking requirements') # TODO: Get rid of!
# check_requirements()
print('write version.py')
write_version_py()
setup(name=DISTNAME,
description=DESCRIPTION,
long_description=LONG_DESCRIPTION,
maintainer=MAINTAINER,
maintainer_email=MAINTAINER_EMAIL,
url=URL,
download_url=URL,
version=VERSION,
packages=find_packages(exclude=['tests', 'doc']), # TODO: necessary?
#include_dirs=[numpy.get_include()], # TODO: Maybe used for sphinx?!
#setup_requires=['numpy>=1.6', 'pytest', 'pytest-runner'],
#tests_require=['pytest', 'pytest-cov', 'pytest-flake8'],
install_requires=install_requires,
extras_require=extras_require,
#cmdclass={'build': build} # TODO: necessary?
)
print('-------------------------------------------------------------------------------\n')
src/empyre/__init__.py 0 → 100644
View file @ 6e3128b8
# -*- coding: utf-8 -*-
# Copyright 2020 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing functionality for visualisation of multidimensional fields."""
from . import fields
from . import io
from . import models
from . import reconstruct
from . import vis
from . import utils
__version__ = "0.3.4"
__git_revision__ = "undefined"
import logging
_log = logging.getLogger(__name__)
_log.info(f'Imported EMPyRe V-{__version__}')
del logging
__all__ = ['fields', 'io', 'models', 'reconstruct', 'vis', 'utils']